diff --git a/cosserat-continuum-twisted-strip.parset b/cosserat-continuum-twisted-strip.parset
index d2cdfd399e11c2e54fc2b7e9e10432035ce508c3..6accef9b1ac413e15fc749e15074a6622ec7b15e 100644
--- a/cosserat-continuum-twisted-strip.parset
+++ b/cosserat-continuum-twisted-strip.parset
@@ -84,6 +84,8 @@ kappa = 1
# Boundary values
#############################################
+problem = twisted-strip
+
### Python predicate specifying all Dirichlet grid vertices
# x is the vertex coordinate
dirichletVerticesPredicate = "x[0] < 0.001 or x[0] > 0.0999"
diff --git a/cosserat-continuum-wong-pellegrino-wrinkling.parset b/cosserat-continuum-wong-pellegrino-wrinkling.parset
index e01e58cf48f2d2999e280265932dd88dc289af22..790a91cd3f649131399a63fbccc924a41ee268c5 100644
--- a/cosserat-continuum-wong-pellegrino-wrinkling.parset
+++ b/cosserat-continuum-wong-pellegrino-wrinkling.parset
@@ -21,7 +21,7 @@ numHomotopySteps = 1
tolerance = 1e-8
# Max number of steps of the trust region solver
-maxTrustRegionSteps = 200
+maxTrustRegionSteps = 500
# Initial trust-region radius
initialTrustRegionRadius = 0.001
@@ -71,7 +71,7 @@ lambda = 2.1796e+09
mu_c = 0
# Length scale parameter
-L_c = 2.5e-5
+L_c = 2.5e-8
# Curvature exponent
q = 2
@@ -92,4 +92,6 @@ problem = wong-pellegrino
dirichletVerticesPredicate = "x[1] < 0.0001 or x[1] > 0.128 - 0.0001"
# Initial deformation
+#startFromFile = true
+initialIterateFilename = cosserat_iterate_2.vtu
initialDeformation = "[x[0] + 0.003*x[1] / 0.128, x[1], 0.002*math.cos(1e4*x[0])]"
diff --git a/cosserat-continuum-wriggers-l-shape.parset b/cosserat-continuum-wriggers-l-shape.parset
index d9be4481e201a77898a80d4aaf1d78e1c60a3dac..73364963271993ea68c7e5b65517b7b4ebd0ea61 100644
--- a/cosserat-continuum-wriggers-l-shape.parset
+++ b/cosserat-continuum-wriggers-l-shape.parset
@@ -17,10 +17,10 @@ numLevels = 1
numHomotopySteps = 1
# Tolerance of the trust region solver
-tolerance = 1e-8
+tolerance = 1e-6
# Max number of steps of the trust region solver
-maxTrustRegionSteps = 200
+maxTrustRegionSteps = 500
# Initial trust-region radius
initialTrustRegionRadius = 0.1
@@ -76,7 +76,7 @@ L_c = 0.6e-3
q = 2
# Shear correction factor
-kappa = 1
+kappa = 0.01
[]
@@ -84,6 +84,8 @@ kappa = 1
# Boundary values
#############################################
+problem = wriggers-l-shape
+
### Python predicate specifying all Dirichlet grid vertices
# x is the vertex coordinate
dirichletVerticesPredicate = "x[0] < 0.001"
@@ -93,7 +95,8 @@ dirichletVerticesPredicate = "x[0] < 0.001"
neumannVerticesPredicate = "x[1] < -0.239"
### Neumann values, if needed
-neumannValues = -5e6 0 0
+#neumannValues = -2e7 0 0
# Initial deformation
-initialDeformation = "[x[0], x[1], -0.1*x[0]*(0.24+x[1])]"
+#initialDeformation = "[x[0], x[1], -0.3*x[0]*(0.24+x[1])]"
+initialDeformation = "[x[0], x[1], 0 if (x[0] < 0.225 or x[1] <-0.015) else 4*(x[0]-0.225)*(x[1]+0.015)]"
diff --git a/dune/gfe/trustregionmmgbasesolver.hh b/dune/gfe/trustregionmmgbasesolver.hh
index 48db1958eb718fb64894f264c2161e7e16e5df01..7b4d9aa29f193af99ca5841e4314b168e0b91cde 100644
--- a/dune/gfe/trustregionmmgbasesolver.hh
+++ b/dune/gfe/trustregionmmgbasesolver.hh
@@ -3,6 +3,7 @@
#ifndef DUNE_GFE_TRUSTREGIONMMGBASESOLVER_HH
#define DUNE_GFE_TRUSTREGIONMMGBASESOLVER_HH
+#include <dune/istl/solver.hh>
#include <dune/istl/umfpack.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
@@ -101,6 +102,10 @@ void TrustRegionMMGBaseSolver<MatrixType, VectorType>::solve()
// Modify Dirichlet rows
///////////////////////////////////////////
+ for (size_t j=0; j<modifiedStiffness.N(); j++)
+ for (int k=0; k<blocksize; k++)
+ modifiedStiffness[j][j][k][k] += 0.0;
+
for (size_t j=0; j<modifiedStiffness.N(); j++)
{
auto cIt = modifiedStiffness[j].begin();
@@ -128,6 +133,7 @@ void TrustRegionMMGBaseSolver<MatrixType, VectorType>::solve()
/////////////////////////////////////////////////////////////////
Dune::InverseOperatorResult statistics;
Dune::UMFPack<MatrixType> solver(modifiedStiffness);
+ solver.setOption(UMFPACK_PRL, 0); // no output at all
solver.apply(*x_, modifiedRhs, statistics);
// Model increase? Then the matrix is not positive definite -- fall back to ipopt solver
@@ -138,17 +144,21 @@ void TrustRegionMMGBaseSolver<MatrixType, VectorType>::solve()
matrix_->umv(*x_, tmp);
double modelDecrease = (modifiedRhs*(*x_)) - 0.5 * ((*x_)*tmp);
- if (std::isnan(modelDecrease) or modelDecrease < 0)
+ if (std::isnan(modelDecrease) or modelDecrease < 0 or true)
{
std::cout << "Model increase: " << -modelDecrease << ", falling back to slower solver" << std::endl;
+ std::cout << "Total VARIABLES: " << x_->size() << " x " << blocksize << " = " << x_->size()*blocksize << std::endl;
+ std::cout << "Dirichlet DOFS: " << this->ignoreNodes_->count() << std::endl;
+
QuadraticIPOptSolver<MatrixType, VectorType> baseSolver;
- baseSolver.verbosity_ = NumProc::QUIET;
+ baseSolver.verbosity_ = NumProc::REDUCED;
baseSolver.tolerance_ = 1e-8;
*x_ = 0;
baseSolver.setProblem(*matrix_, *x_, *rhs_);
baseSolver.obstacles_ = obstacles_;
+ baseSolver.ignoreNodes_ = this->ignoreNodes_;
baseSolver.solve();
}
@@ -173,12 +183,13 @@ void TrustRegionMMGBaseSolver<MatrixType, VectorType>::solve()
std::cout << "Inadmissible step, falling back to slower solver" << std::endl;
QuadraticIPOptSolver<MatrixType, VectorType> baseSolver;
- baseSolver.verbosity_ = NumProc::QUIET;
+ baseSolver.verbosity_ = NumProc::REDUCED;
baseSolver.tolerance_ = 1e-8;
*x_ = 0;
baseSolver.setProblem(*matrix_, *x_, *rhs_);
baseSolver.obstacles_ = obstacles_;
+ baseSolver.ignoreNodes_ = this->ignoreNodes_;
baseSolver.solve();
}
diff --git a/finite-strain-elasticity.cc b/finite-strain-elasticity.cc
index 4f99570349a934e6e544e64e9ffca590f8dc53e6..512449410f1a9a1a83b0883545f953ce40806838 100644
--- a/finite-strain-elasticity.cc
+++ b/finite-strain-elasticity.cc
@@ -134,7 +134,8 @@ int main (int argc, char *argv[]) try
typedef GridType::LeafGridView GridView;
GridView gridView = grid->leafGridView();
- typedef P1NodalBasis<GridView,double> FEBasis;
+// typedef P1NodalBasis<GridView,double> FEBasis;
+ typedef P2NodalBasis<GridView,double> FEBasis;
FEBasis feBasis(gridView);
// /////////////////////////////////////////
@@ -271,14 +272,14 @@ int main (int argc, char *argv[]) try
std::vector<FieldVector<double,3> > pointLoads(x.size());
std::fill(pointLoads.begin(), pointLoads.end(), 0);
- pointLoads[1372] = parameterSet.get<FieldVector<double,3> >("neumannValues");
- pointLoads[1372] *= 0.5;
+// pointLoads[1372] = parameterSet.get<FieldVector<double,3> >("neumannValues");
+// pointLoads[1372] *= 0.5;
// /////////////////////////////////////////////////
// Create a Riemannian trust-region solver
// /////////////////////////////////////////////////
- TrustRegionSolver<GridType,SolutionType> solver;
+ TrustRegionSolver<FEBasis,SolutionType> solver;
solver.setup(*grid,
&assembler,
x,
diff --git a/harmonicmaps.cc b/harmonicmaps.cc
index d37a1c3ae57ce11f118ca4a6ea756cbbd6443121..d8de45aa40532062881d9f32cd25d5301e72ae5d 100644
--- a/harmonicmaps.cc
+++ b/harmonicmaps.cc
@@ -33,6 +33,7 @@
#include <dune/gfe/unitvector.hh>
#include <dune/gfe/realtuple.hh>
#include <dune/gfe/harmonicenergystiffness.hh>
+#include <dune/gfe/chiralskyrmionenergy.hh>
#include <dune/gfe/geodesicfeassembler.hh>
#include <dune/gfe/riemanniantrsolver.hh>
@@ -229,6 +230,7 @@ int main (int argc, char *argv[]) try
typedef P1NodalBasis<typename GridType::LeafGridView,double> FEBasis;
FEBasis feBasis(grid.leafGridView());
+ GFE::ChiralSkyrmionEnergy<GridType::LeafGridView, FEBasis::LocalFiniteElement, double> chiralSkyrmionEnergy;
HarmonicEnergyLocalStiffness<GridType::LeafGridView, FEBasis::LocalFiniteElement, TargetSpace> harmonicEnergyLocalStiffness;
GeodesicFEAssembler<FEBasis,TargetSpace> assembler(grid.leafGridView(),
diff --git a/mixed-cosserat-continuum.cc b/mixed-cosserat-continuum.cc
index 376cb2bef8cb236a9d75b0017b7605b13257ae90..993b0582167014fc2c7b7dc7069f072dbf102372 100644
--- a/mixed-cosserat-continuum.cc
+++ b/mixed-cosserat-continuum.cc
@@ -243,12 +243,14 @@ int main (int argc, char *argv[]) try
typedef GridType::LeafGridView GridView;
GridView gridView = grid->leafGridView();
- typedef P2NodalBasis<GridView,double> DeformationFEBasis;
- typedef P1NodalBasis<GridView,double> OrientationFEBasis;
+ typedef P3NodalBasis<GridView,double> DeformationFEBasis;
+ typedef P2NodalBasis<GridView,double> OrientationFEBasis;
DeformationFEBasis deformationFEBasis(gridView);
OrientationFEBasis orientationFEBasis(gridView);
+ std::cout << "Deformation: " << deformationFEBasis.size() << ", orientation: " << orientationFEBasis.size() << std::endl;
+
// /////////////////////////////////////////
// Read Dirichlet values
// /////////////////////////////////////////
diff --git a/test/Makefile.am b/test/Makefile.am
index adf1c5c5959cbee7a3d0cd307aff88de458104a6..d56b4d3641b2b267fd28d3719de6a298f5a66a18 100644
--- a/test/Makefile.am
+++ b/test/Makefile.am
@@ -19,6 +19,7 @@ check_PROGRAMS = adolctest \
globalgfetestfunctionbasistest \
harmonicenergytest \
interpolationtest \
+ interillustration \
localgeodesicfefunctiontest \
localgeodesicfestiffnesstest \
localgfetestfunctiontest \
@@ -30,11 +31,12 @@ check_PROGRAMS = adolctest \
rotationtest \
svdtest \
targetspacetest \
- true-adolctest
+ true-adolctest \
+ vtkreadertest
adolctest_SOURCES = adolctest.cc
adolctest_LDFLAGS = $(ADOLC_LDFLAGS) $(PYTHON_LDFLAGS) $(IPOPT_LDFLAGS)
-adolctest_LDADD = $(ADOLC_LIBS) $(PYTHON_LIBS) $(IPOPT_LIBS)
+adolctest_LDADD = $(ADOLC_LIBS) $(PYTHON_LIBS) $(IPOPT_LIBS) -lmpfr -lgmpxx -lgmp
true_adolctest_SOURCES = true-adolctest.cc
true_adolctest_LDFLAGS = $(ADOLC_LDFLAGS) $(PYTHON_LDFLAGS) $(IPOPT_LDFLAGS)
@@ -64,6 +66,8 @@ harmonicenergytest_SOURCES = harmonicenergytest.cc
interpolationtest_SOURCES = interpolationtest.cc
+interillustration_SOURCES = interillustration.cc
+
cosseratenergytest_SOURCES = cosseratenergytest.cc
averagedistanceassemblertest_SOURCES = averagedistanceassemblertest.cc
@@ -76,6 +80,10 @@ targetspacetest_SOURCES = targetspacetest.cc
svdtest_SOURCES = svdtest.cc
+vtkreadertest_SOURCES = vtkreadertest.cc
+vtkreadertest_CXXFLAGS = -DHAVE_TINYXML2
+vtkreadertest_LDADD = -ltinyxml2
+
# don't follow the full GNU-standard
# we need automake 1.5
AUTOMAKE_OPTIONS = foreign 1.5
diff --git a/test/adolctest.cc b/test/adolctest.cc
index 2e60666a81a113a8e5fd34ea955690b992221642..bfbcc17f18d8e3208727ab29b5703dbf914290fe 100644
--- a/test/adolctest.cc
+++ b/test/adolctest.cc
@@ -30,6 +30,7 @@ typedef double FDType;
#include <dune/istl/io.hh>
#include <dune/fufem/functionspacebases/p2nodalbasis.hh>
+#include <dune/fufem/functiontools/basisinterpolator.hh>
#include <dune/gfe/rigidbodymotion.hh>
@@ -47,7 +48,19 @@ typedef RigidBodyMotion<double,3> TargetSpace;
using namespace Dune;
+class Identity
+: public VirtualFunction<FieldVector<double,2>, FieldVector<double,3> >
+{
+public:
+ void evaluate(const FieldVector<double,2>& in,
+ FieldVector<double,3>& out) const
+ {
+ out[0] = in[0];
+ out[1] = in[1];
+ out[2] = 0.0;
+ }
+};
/** \brief Assembles energy gradient and Hessian with ADOL-C
*/
@@ -446,6 +459,7 @@ int main (int argc, char *argv[]) try
//////////////////////////////////////////7
// Read initial iterate from file
//////////////////////////////////////////7
+#if 0
Dune::BlockVector<FieldVector<double,7> > xEmbedded(x.size());
std::ifstream file("dangerous_iterate", std::ios::in|std::ios::binary);
@@ -458,6 +472,15 @@ int main (int argc, char *argv[]) try
for (int ii=0; ii<x.size(); ii++)
x[ii] = xEmbedded[ii];
+#else
+ Identity identity;
+
+ std::vector<FieldVector<double,3> > v;
+ Functions::interpolate(feBasis, v, identity);
+
+ for (size_t i=0; i<x.size(); i++)
+ x[i].r = v[i];
+#endif
// ////////////////////////////////////////////////////////////
// Create an assembler for the energy functional
diff --git a/test/targetspacetest.cc b/test/targetspacetest.cc
index 67aa95884d4ae559f1eef8d63d6c3cd0a24dc467..820f3e81c8199abf86be9a924fe6ecbac1bb5af0 100644
--- a/test/targetspacetest.cc
+++ b/test/targetspacetest.cc
@@ -345,7 +345,7 @@ void test()
// Test each element in the list
for (int i=0; i<nTestPoints; i++) {
- testOrthonormalFrame<TargetSpace>(testPoints[i]);
+ //testOrthonormalFrame<TargetSpace>(testPoints[i]);
for (int j=0; j<nTestPoints; j++) {
@@ -354,8 +354,13 @@ void test()
testPointPair[1] = testPoints[j];
if (diameter(testPointPair) > TargetSpace::convexityRadius)
continue;
+
+ TargetSpace p = testPointPair[0];
+ TargetSpace q = testPointPair[1];
+ std::cout << "p: " << testPointPair[0] << ", q: " << testPointPair[1] << std::endl;
+ std::cout << TargetSpace::exp(p, TargetSpace::log(p,q)) << std::endl;
- testDerivativesOfSquaredDistance<TargetSpace>(testPoints[i], testPoints[j]);
+ //testDerivativesOfSquaredDistance<TargetSpace>(testPoints[i], testPoints[j]);
}
@@ -366,18 +371,18 @@ void test()
int main() try
{
- test<RealTuple<double,1> >();
- test<RealTuple<double,3> >();
+// test<RealTuple<double,1> >();
+// test<RealTuple<double,3> >();
test<UnitVector<double,2> >();
test<UnitVector<double,3> >();
test<UnitVector<double,4> >();
- test<Rotation<double,3> >();
-
- test<RigidBodyMotion<double,3> >();
-
- test<HyperbolicHalfspacePoint<double,2> >();
+// test<Rotation<double,3> >();
+//
+// test<RigidBodyMotion<double,3> >();
+//
+// test<HyperbolicHalfspacePoint<double,2> >();
} catch (Exception e) {