diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index e2e965b05c5e080eccab10e7f19aed93aca2eda1..b3c945cd28d141a75cdc5ac92ed1c07a1656c562 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -59,7 +59,7 @@ dune:git clang:
- *before
script: duneci-standard-test
-dune:git dune-parmg dune-vtk dune-curvedgeometry dune-curvedgrid gcc:
+dune:git dune-parmg dune-vtk dune-curvedgeometry dune-curvedgrid dune-foamgrid gcc:
image: registry.dune-project.org/docker/ci/dune:git-debian-10-gcc-8-17
before_script:
- *patch-dune-common
@@ -68,9 +68,10 @@ dune:git dune-parmg dune-vtk dune-curvedgeometry dune-curvedgrid gcc:
- duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/spraetor/dune-vtk.git
- duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/spraetor/dune-curvedgeometry.git
- duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/iwr/dune-curvedgrid.git
+ - duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git
script: duneci-standard-test
-dune:git dune-parmg dune-vtk dune-curvedgeometry dune-curvedgrid clang:
+dune:git dune-parmg dune-vtk dune-curvedgeometry dune-curvedgrid dune-foamgrid clang:
image: registry.dune-project.org/docker/ci/dune:git-ubuntu-20.04-clang-10-20
before_script:
- *patch-dune-common
@@ -79,6 +80,7 @@ dune:git dune-parmg dune-vtk dune-curvedgeometry dune-curvedgrid clang:
- duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/spraetor/dune-vtk.git
- duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/spraetor/dune-curvedgeometry.git
- duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/iwr/dune-curvedgrid.git
+ - duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git
script: duneci-standard-test
# Check for spelling mistakes in text
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 5f2f6dd8e3609ae8df829bd6d8633349c52af710..861e472bc49ad61239d78def7b35a43340bfd581 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,11 @@
# Master
+- Build cosserat-continuum for different combinations of LFE-orders and
+ GFE-orders, the respective program is called
+ cosserat-continuum-Xd-in-Xd-LFE_ORDER-GFE_ORDER
+ TODO: This is now set during compile time, but shall be changed to be
+ set during runtime.
+
- Do not scale the density functions with the thickness to avoid confusion
since some densities need to be scaled and some do not need to be scaled
with the thickness depending on the dimension of the grid, their direction
diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index 2f1eb8b27af6b78dcc48e9a7f4dd14843df4965b..df4749eb9be221a20381c136f34bfa3218e99013 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -21,7 +21,13 @@
#define DONT_USE_CURL
-//#define QUADRATIC_MEMBRANE_ENERGY
+#define CURVATURE_WITH_WRYNESS //only relevant for gridDim == 3
+
+//#define QUADRATIC_2006 //only relevant for gridDim == 3
+// Use the formula for the quadratic membrane energy (first equation of (4.4)) in Neff's paper from 2006: A geometrically exact planar Cosserat shell model with microstructure: Existence of minimizers for zero Cosserat couple modulus
+// instead of the equation (2.27) of 2019: Reï¬ned dimensional reduction for isotropic elastic Cosserat shells with initial curvature
+
+//#define QUADRATIC_MEMBRANE_ENERGY //only relevant for gridDim == 2
/** \brief Get LocalFiniteElements from a localView, for different tree depths of the local view
*
@@ -102,7 +108,7 @@ public:
neumannFunction_(neumannFunction),
volumeLoad_(volumeLoad)
{
- // The shell thickness
+ // The shell thickness // only relevant for dim == 2
thickness_ = parameters.template get<double>("thickness");
// Lame constants
@@ -118,8 +124,13 @@ public:
// Curvature exponent
q_ = parameters.template get<double>("q");
- // Shear correction factor
+ // Shear correction factor // only relevant for dim == 2
kappa_ = parameters.template get<double>("kappa");
+
+ // Curvature parameters
+ b1_ = parameters.template get<double>("b1");
+ b2_ = parameters.template get<double>("b2");
+ b3_ = parameters.template get<double>("b3");
}
/** \brief Assemble the energy for a single element */
@@ -133,6 +144,7 @@ public:
/** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
* the first equation of (4.4) in Neff's paper from 2006: A geometrically exact planar Cosserat shell model with microstructure: Existence of minimizers for zero Cosserat couple modulus
+ * OR: the equation (2.27) of 2019: Reï¬ned dimensional reduction for isotropic elastic Cosserat shells with initial curvature
*/
RT quadraticMembraneEnergy(const Dune::GFE::CosseratStrain<field_type,3,gridDim>& U) const
{
@@ -142,7 +154,11 @@ public:
return mu_ * Dune::GFE::sym(UMinus1).frobenius_norm2()
+ mu_c_ * Dune::GFE::skew(UMinus1).frobenius_norm2()
+#ifdef QUADRATIC_2006
+ (mu_*lambda_)/(2*mu_ + lambda_) * Dune::GFE::traceSquared(UMinus1); // Dune::GFE::traceSquared(UMinus1) = Dune::GFE::traceSquared(Dune::GFE::sym(UMinus1))
+#else
+ + lambda_/2 * Dune::GFE::traceSquared(UMinus1); // Dune::GFE::traceSquared(UMinus1) = Dune::GFE::traceSquared(Dune::GFE::sym(UMinus1))
+#endif
}
/** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
@@ -236,6 +252,36 @@ public:
#endif
}
+ RT curvatureWithWryness(const Dune::FieldMatrix<field_type,dim,dim>& R, const Tensor3<field_type,3,3,gridDim>& DR) const
+ {
+ // construct Wryness tensor \Gamma as in "Refined dimensional reduction for isotropic elastic Cosserat shells with initial curvature"
+ Dune::FieldMatrix<field_type,3,3> dRx1(0); //Derivative of R with respect to x1
+ Dune::FieldMatrix<field_type,3,3> dRx2(0); //Derivative of R with respect to x2
+ Dune::FieldMatrix<field_type,3,3> dRx3(0); //Derivative of R with respect to x3, this of course only exists if gridDim = 3 (then dim = 3 also, because dim >= gridDim, always)
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++) {
+ dRx1[i][j] = DR[i][j][0];
+ dRx2[i][j] = DR[i][j][1];
+ dRx3[i][j] = gridDim == 3 ? DR[i][j][2] : 0;
+ }
+
+ Dune::FieldMatrix<field_type,3,3> wryness(0);
+
+ auto axialVectorx1 = SkewMatrix<field_type,3>(transpose(R)*dRx1).axial();
+ auto axialVectorx2 = SkewMatrix<field_type,3>(transpose(R)*dRx2).axial();
+ auto axialVectorx3 = SkewMatrix<field_type,3>(transpose(R)*dRx3).axial();
+ for (int i=0; i<3; i++) {
+ wryness[i][0] = axialVectorx1[i];
+ wryness[i][1] = axialVectorx2[i];
+ wryness[i][2] = gridDim == 3 ? axialVectorx3[i] : 0;
+ }
+
+ // The choice for the Frobenius norm here is b1=b2=1 and b3 = 1/3
+ return mu_ * L_c_ * L_c_ * (b1_ * Dune::GFE::dev(Dune::GFE::sym(wryness)).frobenius_norm2()
+ + b2_ * Dune::GFE::skew(wryness).frobenius_norm2() + b3_ * Dune::GFE::traceSquared(wryness));
+ }
+
+
RT bendingEnergy(const Dune::FieldMatrix<field_type,dim,dim>& R, const Tensor3<field_type,3,3,gridDim>& DR) const
{
// left-multiply the derivative of the third director (in DR[][2][]) with R^T
@@ -268,6 +314,9 @@ public:
/** \brief Shear correction factor */
double kappa_;
+ /** \brief Curvature parameters */
+ double b1_, b2_, b3_;
+
/** \brief The Neumann boundary */
const BoundaryPatch<GridView>* neumannBoundary_;
@@ -357,7 +406,11 @@ energy(const typename Basis::LocalView& localView,
energy += weight * std::pow(thickness_,3) / 12.0 * bendingEnergy(R,DR);
} else if (gridDim==3) {
energy += weight * quadraticMembraneEnergy(U);
+#ifdef CURVATURE_WITH_WRYNESS
+ energy += weight * curvatureWithWryness(R,DR);
+#else
energy += weight * curvatureEnergy(DR);
+#endif
} else
DUNE_THROW(Dune::NotImplemented, "CosseratEnergyStiffness for 1d grids");
@@ -511,7 +564,11 @@ energy(const typename Basis::LocalView& localView,
energy += weight * std::pow(thickness_,3) / 12.0 * bendingEnergy(R,DR);
} else if (gridDim==3) {
energy += weight * quadraticMembraneEnergy(U);
+#ifdef CURVATURE_WITH_WRYNESS
+ energy += weight * curvatureWithWryness(R,DR);
+#else
energy += weight * curvatureEnergy(DR);
+#endif
} else
DUNE_THROW(Dune::NotImplemented, "CosseratEnergyStiffness for 1d grids");
diff --git a/dune/gfe/nonplanarcosseratshellenergy.hh b/dune/gfe/nonplanarcosseratshellenergy.hh
index cfa82d662c2cfe22508d198cad580e304f978295..4901f8d5f6d1aaf08437c28ce1ca716d36496fac 100644
--- a/dune/gfe/nonplanarcosseratshellenergy.hh
+++ b/dune/gfe/nonplanarcosseratshellenergy.hh
@@ -135,6 +135,7 @@ public:
return mu_ * Dune::GFE::sym(S).frobenius_norm2() + mu_c_ * Dune::GFE::skew(S).frobenius_norm2() + lambda_ * 0.5 * Dune::GFE::traceSquared(S);
}
+ //For b1 = b2 = 1 and b3 = 1/3, this reduces to S.frobenius_norm2()
RT W_curv(const Dune::FieldMatrix<field_type,3,3>& S) const
{
return mu_ * L_c_ * L_c_ * (b1_ * Dune::GFE::dev(Dune::GFE::sym(S)).frobenius_norm2()
diff --git a/dune/gfe/riemannianpnsolver.cc b/dune/gfe/riemannianpnsolver.cc
index a10252d6bbc6e930d7361581eeeabaea106e84af..1085ad501533d53e70d075e3d10a7da3c0fee2e7 100644
--- a/dune/gfe/riemannianpnsolver.cc
+++ b/dune/gfe/riemannianpnsolver.cc
@@ -21,7 +21,8 @@
#include <dune/solvers/norms/twonorm.hh>
#include <dune/solvers/norms/h1seminorm.hh>
-#if HAVE_MPI
+// The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
+#if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3)) && (!defined(WORLD_DIM) or (defined(WORLD_DIM) && WORLD_DIM < 3))
#include <dune/gfe/parallel/matrixcommunicator.hh>
#include <dune/gfe/parallel/vectorcommunicator.hh>
#endif
@@ -46,15 +47,13 @@ setup(const GridType& grid,
instrumented_ = instrumented;
ignoreNodes_ = &dirichletNodes;
+// The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
+#if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3)) && (!defined(WORLD_DIM) or (defined(WORLD_DIM) && WORLD_DIM < 3))
//////////////////////////////////////////////////////////////////
// Create global numbering for matrix and vector transfer
//////////////////////////////////////////////////////////////////
-#if HAVE_MPI
globalMapper_ = std::make_unique<GlobalMapper>(grid_->leafGridView());
-#endif
-
-#if HAVE_MPI
// Transfer all Dirichlet data to the master processor
VectorCommunicator<GlobalMapper, typename GridType::LeafGridView::CollectiveCommunication, Dune::BitSetVector<blocksize> > vectorComm(*globalMapper_,
grid_->leafGridView().comm(),
@@ -77,7 +76,8 @@ setup(const GridType& grid,
operatorAssembler.assembleBulk(Dune::Fufem::istlMatrixBackend(localA), laplaceStiffness);
-#if HAVE_MPI
+// The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
+#if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3)) && (!defined(WORLD_DIM) or (defined(WORLD_DIM) && WORLD_DIM < 3))
LocalMapper localMapper = MapperFactory<Basis>::createLocalMapper(grid_->leafGridView());
MatrixCommunicator<GlobalMapper,
@@ -115,7 +115,8 @@ setup(const GridType& grid,
operatorAssembler.assembleBulk(Dune::Fufem::istlMatrixBackend(localMassMatrix), massStiffness);
-#if HAVE_MPI
+// The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
+#if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3)) && (!defined(WORLD_DIM) or (defined(WORLD_DIM) && WORLD_DIM < 3))
auto massMatrix = std::make_shared<ScalarMatrixType>(matrixComm.reduceAdd(localMassMatrix));
#else
auto massMatrix = std::make_shared<ScalarMatrixType>(localMassMatrix);
@@ -170,7 +171,8 @@ void RiemannianProximalNewtonSolver<Basis,TargetSpace,Assembler>::solve()
bool recomputeGradientHessian = true;
CorrectionType rhs, rhs_global;
MatrixType stiffnessMatrix;
-#if HAVE_MPI
+// The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
+#if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3)) && (!defined(WORLD_DIM) or (defined(WORLD_DIM) && WORLD_DIM < 3))
VectorCommunicator<GlobalMapper, typename GridType::LeafGridView::CollectiveCommunication, CorrectionType> vectorComm(*globalMapper_,
grid_->leafGridView().comm(),
0);
@@ -216,7 +218,8 @@ void RiemannianProximalNewtonSolver<Basis,TargetSpace,Assembler>::solve()
rhs *= -1; // The right hand side is the _negative_ gradient
// Transfer vector data
-#if HAVE_MPI
+// The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
+#if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3)) && (!defined(WORLD_DIM) or (defined(WORLD_DIM) && WORLD_DIM < 3))
rhs_global = vectorComm.reduceAdd(rhs);
#else
rhs_global = rhs;
@@ -234,8 +237,8 @@ void RiemannianProximalNewtonSolver<Basis,TargetSpace,Assembler>::solve()
std::cout << "Overall assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
totalAssemblyTime += gradientTimer.elapsed();
- // Transfer matrix data
-#if HAVE_MPI
+// The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
+#if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3)) && (!defined(WORLD_DIM) or (defined(WORLD_DIM) && WORLD_DIM < 3))
stiffnessMatrix = matrixComm.reduceAdd(*hessianMatrix_);
#else
stiffnessMatrix = *hessianMatrix_;
@@ -253,7 +256,7 @@ void RiemannianProximalNewtonSolver<Basis,TargetSpace,Assembler>::solve()
// Add the regularization - Identity Matrix for now
for (std::size_t i=0; i<stiffnessMatrix.N(); i++)
for(int j=0; j<blocksize; j++)
- stiffnessMatrix[i][i][j][j] += regularization;
+ stiffnessMatrix[i][i][j][j] += regularization*scaling_[j];
innerSolver_->setProblem(stiffnessMatrix,corr_global,rhs_global);
innerSolver_->preprocess();
@@ -280,7 +283,8 @@ void RiemannianProximalNewtonSolver<Basis,TargetSpace,Assembler>::solve()
if (grid_->comm().size()>1 and rank==0)
std::cout << "Transfer solution back to root process ..." << std::endl;
-#if HAVE_MPI
+// The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
+#if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3)) && (!defined(WORLD_DIM) or (defined(WORLD_DIM) && WORLD_DIM < 3))
solved = grid_->comm().min(solved);
if (solved) {
corr = vectorComm.scatter(corr_global);
diff --git a/dune/gfe/riemannianpnsolver.hh b/dune/gfe/riemannianpnsolver.hh
index b3584fdf9b8553e0311a671d8e3b14732fd3e5f8..7ea10fe6f4737bf6005560134cf52e3f281338d4 100644
--- a/dune/gfe/riemannianpnsolver.hh
+++ b/dune/gfe/riemannianpnsolver.hh
@@ -69,7 +69,9 @@ public:
RiemannianProximalNewtonSolver()
: IterativeSolver<std::vector<TargetSpace>, Dune::BitSetVector<blocksize> >(0,100,NumProc::FULL),
hessianMatrix_(nullptr), h1SemiNorm_(NULL)
- {}
+ {
+ std::fill(scaling_.begin(), scaling_.end(), 1.0);
+ }
/** \brief Set up the solver using a choldmod or umfpack solver as the inner solver */
void setup(const GridType& grid,
@@ -81,6 +83,11 @@ public:
double initialRegularization,
bool instrumented);
+ void setScaling(const Dune::FieldVector<double,blocksize>& scaling)
+ {
+ scaling_ = scaling;
+ }
+
void setIgnoreNodes(const Dune::BitSetVector<blocksize>& ignoreNodes)
{
ignoreNodes_ = &ignoreNodes;
@@ -118,6 +125,9 @@ protected:
double initialRegularization_;
double tolerance_;
+ /** \brief Regularization scaling */
+ Dune::FieldVector<double,blocksize> scaling_;
+
/** \brief Maximum number of proximal-newton steps */
std::size_t maxProximalNewtonSteps_;
diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index 87ccf8f51bffaafea22574e406ca68a64d672181..d5bdce74ab503275b860496b21f22fe1c7a8746f 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -22,7 +22,8 @@
#include <dune/solvers/norms/twonorm.hh>
#include <dune/solvers/norms/h1seminorm.hh>
-#if HAVE_MPI
+// The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
+#if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3))
#include <dune/gfe/parallel/matrixcommunicator.hh>
#include <dune/gfe/parallel/vectorcommunicator.hh>
#endif
@@ -60,11 +61,12 @@ setup(const GridType& grid,
int numLevels = grid_->maxLevel()+1;
+// The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
+#if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3)) && (!defined(WORLD_DIM) or (defined(WORLD_DIM) && WORLD_DIM < 3))
//////////////////////////////////////////////////////////////////
// Create global numbering for matrix and vector transfer
//////////////////////////////////////////////////////////////////
-#if HAVE_MPI
globalMapper_ = std::make_unique<GlobalMapper>(grid_->leafGridView());
#endif
@@ -89,7 +91,8 @@ setup(const GridType& grid,
baseNorm,
Solver::QUIET);
#endif
-#if HAVE_MPI
+// The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
+#if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3)) && (!defined(WORLD_DIM) or (defined(WORLD_DIM) && WORLD_DIM < 3))
// Transfer all Dirichlet data to the master processor
VectorCommunicator<GlobalMapper, typename GridType::LeafGridView::CollectiveCommunication, Dune::BitSetVector<blocksize> > vectorComm(*globalMapper_,
grid_->leafGridView().comm(),
@@ -124,7 +127,8 @@ setup(const GridType& grid,
operatorAssembler.assembleBulk(Dune::Fufem::istlMatrixBackend(localA), laplaceStiffness);
-#if HAVE_MPI
+// The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
+#if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3)) && (!defined(WORLD_DIM) or (defined(WORLD_DIM) && WORLD_DIM < 3))
LocalMapper localMapper = MapperFactory<Basis>::createLocalMapper(grid_->leafGridView());
MatrixCommunicator<GlobalMapper,
@@ -156,7 +160,8 @@ setup(const GridType& grid,
operatorAssembler.assembleBulk(Dune::Fufem::istlMatrixBackend(localMassMatrix), massStiffness);
-#if HAVE_MPI
+// The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
+#if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3)) && (!defined(WORLD_DIM) or (defined(WORLD_DIM) && WORLD_DIM < 3))
auto massMatrix = std::make_shared<ScalarMatrixType>(matrixComm.reduceAdd(localMassMatrix));
#else
auto massMatrix = std::make_shared<ScalarMatrixType>(localMassMatrix);
@@ -205,7 +210,9 @@ setup(const GridType& grid,
TransferOperatorType pkToP1TransferMatrix;
assembleGlobalBasisTransferMatrix(pkToP1TransferMatrix,p1Basis,basis);
-#if HAVE_MPI
+
+// The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
+#if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3)) && (!defined(WORLD_DIM) or (defined(WORLD_DIM) && WORLD_DIM < 3))
// If we are on more than 1 processors, join all local transfer matrices on rank 0,
// and construct a single global transfer operator there.
typedef Dune::GlobalP1Mapper<Dune::Functions::LagrangeBasis<typename Basis::GridView,1>> GlobalLeafP1Mapper;
@@ -237,7 +244,8 @@ setup(const GridType& grid,
// Construct the local multigrid transfer matrix
auto newTransferOp = std::make_unique<TruncatedCompressedMGTransfer<CorrectionType>>();
newTransferOp->setup(*grid_,i,i+1);
-#if HAVE_MPI
+// The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
+#if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3)) && (!defined(WORLD_DIM) or (defined(WORLD_DIM) && WORLD_DIM < 3))
// If we are on more than 1 processors, join all local transfer matrices on rank 0,
// and construct a single global transfer operator there.
typedef Dune::Functions::LagrangeBasis<typename GridType::LevelGridView, 1> FEBasis;
@@ -250,7 +258,8 @@ setup(const GridType& grid,
LevelLocalMapper coarseLevelLocalMapper(grid_->levelGridView(i), Dune::mcmgVertexLayout());
#endif
typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType;
-#if HAVE_MPI
+// The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
+#if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3)) && (!defined(WORLD_DIM) or (defined(WORLD_DIM) && WORLD_DIM < 3))
MatrixCommunicator<GlobalLevelP1Mapper,
typename GridType::LevelGridView,
typename GridType::LevelGridView,
@@ -274,7 +283,8 @@ setup(const GridType& grid,
if (rank==0)
{
-#if HAVE_MPI
+// The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
+#if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3)) && (!defined(WORLD_DIM) or (defined(WORLD_DIM) && WORLD_DIM < 3))
hasObstacle_.resize(globalMapper_->size(), true);
#else
hasObstacle_.resize(basis.size(), true);
@@ -298,7 +308,8 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::solve()
mgStep = dynamic_cast<MonotoneMGStep<MatrixType,CorrectionType>*>(&loopSolver->getIterationStep());
}
-#if HAVE_MPI
+// The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
+#if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3)) && (!defined(WORLD_DIM) or (defined(WORLD_DIM) && WORLD_DIM < 3))
MaxNormTrustRegion<blocksize> trustRegion(globalMapper_->size(), initialTrustRegionRadius_);
#else
const Basis& basis = assembler_->getBasis();
@@ -328,7 +339,7 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::solve()
Dune::Timer energyTimer;
double oldEnergy = assembler_->computeEnergy(x_);
if (this->verbosity_ == Solver::FULL)
- std::cout << "Energy computation took " << energyTimer.elapsed() << " sec." << std::endl;
+ std::cout << "Energy computation took " << energyTimer.elapsed() << " sec." << ", final energy: " << oldEnergy << std::endl;
oldEnergy = grid_->comm().sum(oldEnergy);
@@ -337,7 +348,8 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::solve()
CorrectionType rhs;
MatrixType stiffnessMatrix;
CorrectionType rhs_global;
-#if HAVE_MPI
+// The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
+#if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3)) && (!defined(WORLD_DIM) or (defined(WORLD_DIM) && WORLD_DIM < 3))
VectorCommunicator<GlobalMapper, typename GridType::LeafGridView::CollectiveCommunication, CorrectionType> vectorComm(*globalMapper_,
grid_->leafGridView().comm(),
0);
@@ -388,7 +400,8 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::solve()
rhs *= -1; // The right hand side is the _negative_ gradient
// Transfer vector data
-#if HAVE_MPI
+// The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
+#if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3)) && (!defined(WORLD_DIM) or (defined(WORLD_DIM) && WORLD_DIM < 3))
rhs_global = vectorComm.reduceAdd(rhs);
#else
rhs_global = rhs;
@@ -407,8 +420,8 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::solve()
std::cout << "Overall assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
totalAssemblyTime += gradientTimer.elapsed();
- // Transfer matrix data
-#if HAVE_MPI
+// The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
+#if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3)) && (!defined(WORLD_DIM) or (defined(WORLD_DIM) && WORLD_DIM < 3))
stiffnessMatrix = matrixComm.reduceAdd(*hessianMatrix_);
#else
stiffnessMatrix = *hessianMatrix_;
@@ -456,7 +469,8 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::solve()
if (grid_->comm().size()>1 and rank==0)
std::cout << "Transfer solution back to root process ..." << std::endl;
-#if HAVE_MPI
+// The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
+#if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3)) && (!defined(WORLD_DIM) or (defined(WORLD_DIM) && WORLD_DIM < 3))
solved = grid_->comm().min(solved);
if (solved) {
corr = vectorComm.scatter(corr_global);
@@ -546,8 +560,8 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::solve()
std::cout << i+1 << " trust-region steps were taken, the maximum was reached." << std::endl << "Total solver time: " << totalSolverTime << " sec., total assembly time: " << totalAssemblyTime << " sec." << std::endl;
if (solved) {
- if (this->verbosity_ == NumProc::FULL)
- std::cout << "Infinity norm of the correction: " << corr.infinity_norm() << std::endl;
+ if (this->verbosity_ == NumProc::FULL and rank==0)
+ std::cout << "Infinity norm of the correction: " << corrGlobalInfinityNorm << std::endl;
if (corrGlobalInfinityNorm < this->tolerance_ && corrGlobalInfinityNorm < trustRegion.radius()*smallestScalingParameter) {
if (this->verbosity_ == NumProc::FULL and rank==0)
diff --git a/problems/cosserat-continuum-cantilever.parset b/problems/cosserat-continuum-cantilever.parset
index 9440301af2c2e3a244b9799310e33768d16e1212..bfa21a8436ebefa2092dc2b059faac5f2d8323fd 100644
--- a/problems/cosserat-continuum-cantilever.parset
+++ b/problems/cosserat-continuum-cantilever.parset
@@ -26,7 +26,7 @@ tolerance = 1e-3
# Max number of steps of the trust region solver
maxSolverSteps = 200
-trustRegionScaling = 1 1 1 0.01 0.01 0.01
+solverScaling = 1 1 1 0.01 0.01 0.01
# Initial trust-region radius
initialTrustRegionRadius = 3.125
diff --git a/problems/cosserat-continuum-twisted-strip.parset b/problems/cosserat-continuum-twisted-strip.parset
index 04c738f4608aa5e30d30c5af4c6a7add32cd6fbe..a10d3197bfb8c791a656639c061b892cb5c48282 100644
--- a/problems/cosserat-continuum-twisted-strip.parset
+++ b/problems/cosserat-continuum-twisted-strip.parset
@@ -23,7 +23,7 @@ tolerance = 1e-8
# Max number of steps of the trust region solver
maxTrustRegionSteps = 200
-trustRegionScaling = 1 1 1 1 1 1
+solverScaling = 1 1 1 1 1 1
# Initial trust-region radius
initialTrustRegionRadius = 0.1
diff --git a/problems/cosserat-continuum-wong-pellegrino-wrinkling.parset b/problems/cosserat-continuum-wong-pellegrino-wrinkling.parset
index 790a91cd3f649131399a63fbccc924a41ee268c5..4af4a38d59003379a994bb2db4d246ac43351ceb 100644
--- a/problems/cosserat-continuum-wong-pellegrino-wrinkling.parset
+++ b/problems/cosserat-continuum-wong-pellegrino-wrinkling.parset
@@ -23,6 +23,8 @@ tolerance = 1e-8
# Max number of steps of the trust region solver
maxTrustRegionSteps = 500
+solverScaling = 1 1 1 0.01 0.01 0.01
+
# Initial trust-region radius
initialTrustRegionRadius = 0.001
diff --git a/problems/cosserat-continuum-wriggers-l-shape.parset b/problems/cosserat-continuum-wriggers-l-shape.parset
index 12c715bd9e9618fa4ce069b8409584294f0f5660..745f333600956d4c13e55af38ef37db5d6e7d74e 100644
--- a/problems/cosserat-continuum-wriggers-l-shape.parset
+++ b/problems/cosserat-continuum-wriggers-l-shape.parset
@@ -22,7 +22,7 @@ tolerance = 1e-8
# Max number of steps of the trust region solver
maxSolverSteps = 1000
-trustRegionScaling = 1 1 1 0.01 0.01 0.01
+solverScaling = 1 1 1 0.01 0.01 0.01
# Initial trust-region radius
initialTrustRegionRadius = 1
@@ -105,4 +105,4 @@ neumannValues = 0.09 0 0
startFromFile = yes
initialIterateGridFilename = wriggers-L-shape_99_mm.msh
-initialIterateFilename = initial-wriggers-l-shape.vtu
\ No newline at end of file
+initialIterateFilename = initial-wriggers-l-shape.vtu
diff --git a/problems/film-on-substrate.parset b/problems/film-on-substrate.parset
index e6cf6ce9164fe7216fc7b8392d764965e02aa2c0..1ce5db4b6676909c3c747f427fddc174dccbc3c7 100644
--- a/problems/film-on-substrate.parset
+++ b/problems/film-on-substrate.parset
@@ -87,7 +87,7 @@ initialRegularization = 1000000
# Solver parameters specific for trust-region solver using multigrid solver
#############################################
-trustRegionScaling = 1 1 1 0.01 0.01 0.01
+solverScaling = 1 1 1 0.01 0.01 0.01
# Initial trust-region radius
initialTrustRegionRadius = 3.125
diff --git a/problems/gradient-flow-curve-shortening.parset b/problems/gradient-flow-curve-shortening.parset
index cf5f31f492d398fe59d167365f950dcb87647e8e..857ec53c5a07417d907c7f09df73b059a5582233 100644
--- a/problems/gradient-flow-curve-shortening.parset
+++ b/problems/gradient-flow-curve-shortening.parset
@@ -20,6 +20,8 @@ tolerance = 1e-8
# Max number of steps of the trust region solver
maxTrustRegionSteps = 100
+solverScaling = 1 1 1 0.01 0.01 0.01
+
# Initial trust-region radius
initialTrustRegionRadius = 0.25
diff --git a/problems/harmonicmaps-skyrmions-hexagon.parset b/problems/harmonicmaps-skyrmions-hexagon.parset
index 1394c0ecd0a14308f4abc54a8fac66e40499044e..4cd7c41b71d19b0e8e76f1183cccf11ec26312e6 100644
--- a/problems/harmonicmaps-skyrmions-hexagon.parset
+++ b/problems/harmonicmaps-skyrmions-hexagon.parset
@@ -21,6 +21,8 @@ tolerance = 1e-12
# Max number of steps of the trust region solver
maxTrustRegionSteps = 100
+solverScaling = 1 1 1 0.01 0.01 0.01
+
# Initial trust-region radius
initialTrustRegionRadius = 1
diff --git a/problems/simofox-cantilever.parset b/problems/simofox-cantilever.parset
index 1eab7e5b40a79956681dc71b7bf0be101e824f3d..4351b50e165c7a9356728eaf4f2049dd80573f8c 100644
--- a/problems/simofox-cantilever.parset
+++ b/problems/simofox-cantilever.parset
@@ -26,7 +26,7 @@ tolerance = 1e-6
# Max number of steps of the trust region solver
maxTrustRegionSteps = 200
-trustRegionScaling = 1 1 1 0.01 0.01
+solverScaling = 1 1 1 0.01 0.01
# Initial trust-region radius
initialTrustRegionRadius = 10
diff --git a/problems/staticrod.parset b/problems/staticrod.parset
index 2d001ecad0fbe15fff3188010c99d570ae8c0909..7f8e5fe57d9b1415fca73ae3aca3effedaf644d2 100644
--- a/problems/staticrod.parset
+++ b/problems/staticrod.parset
@@ -21,6 +21,8 @@ tolerance = 1e-6
# Max number of steps of the trust region solver
maxTrustRegionSteps = 100
+solverScaling = 1 1 1 0.01 0.01 0.01
+
# Initial trust-region radius
initialTrustRegionRadius = 1
diff --git a/problems/zero-volume-load.py b/problems/zero-volume-load.py
new file mode 100644
index 0000000000000000000000000000000000000000..eea7fcfba90db66c4ac5b0b8764878d19d7dd845
--- /dev/null
+++ b/problems/zero-volume-load.py
@@ -0,0 +1,8 @@
+import math
+
+class VolumeLoad:
+ def __init__(self, homotopyParameter):
+ self.homotopyParameter = homotopyParameter
+
+ def volumeLoad(self, x):
+ return [0,0,0]
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 4962aa5d98b6055ca4856b75545ed8a336fda077..ba44aef10080369b9acf10952b3b4bad5d3d1e5a 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,5 +1,4 @@
set(programs compute-disc-error
- cosserat-continuum
gradient-flow
harmonicmaps
simofoxshell)
@@ -10,11 +9,44 @@ endif()
if (${DUNE_ELASTICITY_VERSION} VERSION_GREATER_EQUAL 2.8)
set(programs film-on-substrate-stress-plot ${programs})
endif()
+
+foreach(_program ${programs})
+ add_executable(${_program} ${_program}.cc)
+endforeach()
+
if (dune-parmg_FOUND AND dune-curvedgeometry_FOUND AND ${DUNE_ELASTICITY_VERSION} VERSION_GREATER_EQUAL 2.8)
+ add_executable("film-on-substrate" film-on-substrate.cc)
+ set_property(TARGET "film-on-substrate" APPEND PROPERTY COMPILE_DEFINITIONS "WORLD_DIM=3; LFE_ORDER=2; GFE_ORDER=2")
set(programs film-on-substrate ${programs})
endif()
+
+add_executable("cosserat-continuum-2d-in-2d" cosserat-continuum.cc)
+set_property(TARGET "cosserat-continuum-2d-in-2d" APPEND PROPERTY COMPILE_DEFINITIONS "GRID_DIM=2; WORLD_DIM=2")
+set(programs cosserat-continuum-2d-in-2d ${programs})
+
+if (dune-foamgrid_FOUND)
+ add_executable("cosserat-continuum-2d-in-3d-2-2" cosserat-continuum.cc)
+ set_property(TARGET "cosserat-continuum-2d-in-3d-2-2" APPEND PROPERTY COMPILE_DEFINITIONS "GRID_DIM=2; WORLD_DIM=3; LFE_ORDER=2; GFE_ORDER=2")
+ set(programs cosserat-continuum-2d-in-3d-2-2 ${programs})
+
+ add_executable("cosserat-continuum-2d-in-3d-2-1" cosserat-continuum.cc)
+ set_property(TARGET "cosserat-continuum-2d-in-3d-2-1" APPEND PROPERTY COMPILE_DEFINITIONS "GRID_DIM=2; WORLD_DIM=3; LFE_ORDER=2; GFE_ORDER=1")
+ set(programs cosserat-continuum-2d-in-3d-2-1 ${programs})
+
+ add_executable("cosserat-continuum-2d-in-3d-1-1" cosserat-continuum.cc)
+ set_property(TARGET "cosserat-continuum-2d-in-3d-1-1" APPEND PROPERTY COMPILE_DEFINITIONS "GRID_DIM=2; WORLD_DIM=3; LFE_ORDER=1; GFE_ORDER=1")
+ set(programs cosserat-continuum-2d-in-3d-1-1 ${programs})
+endif()
+
+add_executable("cosserat-continuum-3d-in-3d-2-2" cosserat-continuum.cc)
+set_property(TARGET "cosserat-continuum-3d-in-3d-2-2" APPEND PROPERTY COMPILE_DEFINITIONS "GRID_DIM=3; WORLD_DIM=3; LFE_ORDER=2; GFE_ORDER=2")
+set(programs cosserat-continuum-3d-in-3d-2-2 ${programs})
+
+add_executable("cosserat-continuum-3d-in-3d-2-1" cosserat-continuum.cc)
+set_property(TARGET "cosserat-continuum-3d-in-3d-2-1" APPEND PROPERTY COMPILE_DEFINITIONS "GRID_DIM=3; WORLD_DIM=3; LFE_ORDER=2; GFE_ORDER=1")
+set(programs cosserat-continuum-3d-in-3d-2-1 ${programs})
+
foreach(_program ${programs})
- add_executable(${_program} ${_program}.cc)
target_link_dune_default_libraries(${_program})
add_dune_pythonlibs_flags(${_program})
if (${DUNE_COMMON_VERSION} VERSION_GREATER_EQUAL 2.8)
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index 17333187ba8891fcda1b2ea8faff965ec0f374bd..8de0a6ce0c8fdcbfe1e71c7c1c62ea14da6cf9af 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -1,4 +1,19 @@
-#define MIXED_SPACE 0
+#ifndef LFE_ORDER
+ #define LFE_ORDER 2
+#endif
+
+#ifndef GFE_ORDER
+ #define GFE_ORDER 1
+#endif
+
+#ifndef MIXED_SPACE
+ #if LFE_ORDER != GFE_ORDER
+ #define MIXED_SPACE 1
+ #else
+ #define MIXED_SPACE 0
+ #endif
+#endif
+
//#define PROJECTED_INTERPOLATION
#include <config.h>
@@ -66,14 +81,14 @@
// grid dimension
-const int dim = 2;
-const int dimworld = 2;
+const int dim = GRID_DIM;
+const int dimworld = WORLD_DIM;
// Order of the approximation space for the displacement
-const int displacementOrder = 2;
+const int displacementOrder = LFE_ORDER;
// Order of the approximation space for the microrotations
-const int rotationOrder = 2;
+const int rotationOrder = GFE_ORDER;
#if !MIXED_SPACE
static_assert(displacementOrder==rotationOrder, "displacement and rotation order do not match!");
@@ -106,8 +121,7 @@ int main (int argc, char *argv[]) try
//feenableexcept(FE_INVALID);
Python::runStream()
<< std::endl << "import sys"
- << std::endl << "import os"
- << std::endl << "sys.path.append(os.getcwd() + '/../../problems/')"
+ << std::endl << "sys.path.append('" << argv[1] << "')"
<< std::endl;
using namespace TypeTree::Indices;
@@ -116,10 +130,10 @@ int main (int argc, char *argv[]) try
// parse data file
ParameterTree parameterSet;
- if (argc < 2)
- DUNE_THROW(Exception, "Usage: ./cosserat-continuum <parameter file>");
+ if (argc < 3)
+ DUNE_THROW(Exception, "Usage: ./cosserat-continuum <python path> <parameter file>");
- ParameterTreeParser::readINITree(argv[1], parameterSet);
+ ParameterTreeParser::readINITree(argv[2], parameterSet);
ParameterTreeParser::readOptions(argc, argv, parameterSet);
@@ -139,7 +153,7 @@ int main (int argc, char *argv[]) try
const double baseTolerance = parameterSet.get<double>("baseTolerance");
const bool instrumented = parameterSet.get<bool>("instrumented");
const bool adolcScalarMode = parameterSet.get<bool>("adolcScalarMode", false);
- std::string resultPath = parameterSet.get("resultPath", "");
+ const std::string resultPath = parameterSet.get("resultPath", "");
// ///////////////////////////////////////
// Create the grid
@@ -167,7 +181,7 @@ int main (int argc, char *argv[]) try
grid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);
} else {
- std::string path = parameterSet.get<std::string>("path");
+ std::string path = parameterSet.get<std::string>("path", "");
std::string gridFile = parameterSet.get<std::string>("gridFile");
// Guess the grid file format by looking at the file name suffix
@@ -225,26 +239,45 @@ int main (int argc, char *argv[]) try
)
));
+ BlockVector<FieldVector<double,dimworld> > identityDeformation(compositeBasis.size({0}));
+ auto identityDeformationBasis = makeBasis(
+ gridView,
+ power<dimworld>(
+ lagrange<displacementOrder>()
+ ));
+ Dune::Functions::interpolate(identityDeformationBasis, identityDeformation, [&](FieldVector<double,dimworld> x){ return x;});
+
+
+ BlockVector<FieldVector<double,dimworld> > identityRotation(compositeBasis.size({1}));
+ auto identityRotationBasis = makeBasis(
+ gridView,
+ power<dimworld>(
+ lagrange<rotationOrder>()
+ ));
+ Dune::Functions::interpolate(identityRotationBasis, identityRotation, [&](FieldVector<double,dimworld> x){ return x;});
+
typedef Dune::Functions::LagrangeBasis<GridView,displacementOrder> DeformationFEBasis;
typedef Dune::Functions::LagrangeBasis<GridView,rotationOrder> OrientationFEBasis;
DeformationFEBasis deformationFEBasis(gridView);
OrientationFEBasis orientationFEBasis(gridView);
-
// /////////////////////////////////////////
// Read Dirichlet values
// /////////////////////////////////////////
- BitSetVector<1> dirichletVertices(gridView.size(dim), false);
BitSetVector<1> neumannVertices(gridView.size(dim), false);
+ BitSetVector<3> deformationDirichletDofs(deformationFEBasis.size(), false);
+ BitSetVector<3> orientationDirichletDofs(orientationFEBasis.size(), false);
const GridView::IndexSet& indexSet = gridView.indexSet();
- // Make Python function that computes which vertices are on the Dirichlet boundary,
- // based on the vertex positions.
+ // Make Python function that computes which vertices are on the Dirichlet boundary, based on the vertex positions.
std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletVerticesPredicate") + std::string(")");
- auto pythonDirichletVertices = Python::make_function<bool>(Python::evaluate(lambda));
+ auto pythonDirichletVertices = Python::make_function<FieldVector<bool,3>>(Python::evaluate(lambda));
+
+ lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletRotationVerticesPredicate") + std::string(")");
+ auto pythonOrientationDirichletVertices = Python::make_function<bool>(Python::evaluate(lambda));
// Same for the Neumann boundary
lambda = std::string("lambda x: (") + parameterSet.get<std::string>("neumannVerticesPredicate", "0") + std::string(")");
@@ -252,40 +285,34 @@ int main (int argc, char *argv[]) try
for (auto&& vertex : vertices(gridView))
{
- bool isDirichlet = pythonDirichletVertices(vertex.geometry().corner(0));
- dirichletVertices[indexSet.index(vertex)] = isDirichlet;
-
bool isNeumann = pythonNeumannVertices(vertex.geometry().corner(0));
neumannVertices[indexSet.index(vertex)] = isNeumann;
}
- BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
BoundaryPatch<GridView> neumannBoundary(gridView, neumannVertices);
- std::cout << "On rank " << mpiHelper.rank() << ": Dirichlet boundary has " << dirichletBoundary.numFaces()
- << " faces and " << dirichletVertices.count() << " nodes.\n";
- std::cout << "On rank " << mpiHelper.rank() << ": Neumann boundary has " << neumannBoundary.numFaces() << " faces.\n";
+ std::cout << "On rank " << mpiHelper.rank() << ": Neumann boundary has " << neumannBoundary.numFaces()
+ << " faces and " << neumannVertices.count() << " degrees of freedom.\n";
- BitSetVector<1> deformationDirichletNodes(deformationFEBasis.size(), false);
- constructBoundaryDofs(dirichletBoundary,deformationFEBasis,deformationDirichletNodes);
-
BitSetVector<1> neumannNodes(deformationFEBasis.size(), false);
constructBoundaryDofs(neumannBoundary,deformationFEBasis,neumannNodes);
- BitSetVector<3> deformationDirichletDofs(deformationFEBasis.size(), false);
- for (size_t i=0; i<deformationFEBasis.size(); i++)
- if (deformationDirichletNodes[i][0])
- for (int j=0; j<3; j++)
- deformationDirichletDofs[i][j] = true;
+ for (size_t i=0; i<deformationFEBasis.size(); i++) {
+ FieldVector<bool,3> isDirichlet;
+ isDirichlet = pythonDirichletVertices(identityDeformation[i]);
+ for (size_t j=0; j<3; j++)
+ deformationDirichletDofs[i][j] = isDirichlet[j];
+ }
- BitSetVector<1> orientationDirichletNodes(orientationFEBasis.size(), false);
- constructBoundaryDofs(dirichletBoundary,orientationFEBasis,orientationDirichletNodes);
- BitSetVector<3> orientationDirichletDofs(orientationFEBasis.size(), false);
- for (size_t i=0; i<orientationFEBasis.size(); i++)
- if (orientationDirichletNodes[i][0])
- for (int j=0; j<3; j++)
- orientationDirichletDofs[i][j] = true;
+ for (size_t i=0; i<orientationFEBasis.size(); i++) {
+ bool isDirichletOrientation = pythonOrientationDirichletVertices(identityRotation[i]);
+ for (size_t j=0; j<3; j++)
+ orientationDirichletDofs[i][j] = isDirichletOrientation;
+ }
+
+ std::cout << "On rank " << mpiHelper.rank() << ": Dirichlet boundary has " << deformationDirichletDofs.count() << " degrees of freedom.\n";
+ std::cout << "On rank " << mpiHelper.rank() << ": Dirichlet boundary (orientation) has " << orientationDirichletDofs.count() << " degrees of freedom.\n";
// //////////////////////////
// Initial iterate
@@ -432,14 +459,18 @@ int main (int argc, char *argv[]) try
auto orientationDirichletValues = Python::make_function<FieldMatrix<double,3,3> > (dirichletValuesPythonObject.get("orientation"));
BlockVector<FieldVector<double,3> > ddV;
- Dune::Functions::interpolate(deformationPowerBasis, ddV, deformationDirichletValues, deformationDirichletDofs);
+ Dune::Functions::interpolate(deformationPowerBasis, ddV, deformationDirichletValues);
BlockVector<FieldMatrix<double,3,3> > dOV;
Dune::Functions::interpolate(orientationPowerBasis, dOV, orientationDirichletValues);
for (int i = 0; i < compositeBasis.size({0}); i++) {
- if (deformationDirichletDofs[i][0])
- x[_0][i] = ddV[i];
+ FieldVector<double,3> x0i({x[_0][i][0],x[_0][i][1],x[_0][i][2]});
+ for (int j=0; j<3; j++) {
+ if (deformationDirichletDofs[i][j])
+ x0i[j] = ddV[i][j];
+ }
+ x[_0][i] = x0i;
}
for (int i = 0; i < compositeBasis.size({1}); i++)
if (orientationDirichletDofs[i][0])
@@ -477,7 +508,7 @@ int main (int argc, char *argv[]) try
baseTolerance,
instrumented);
- solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(x);
solver.solve();
@@ -516,7 +547,7 @@ int main (int argc, char *argv[]) try
baseTolerance,
instrumented);
- solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(xTargetSpace);
solver.solve();
xTargetSpace = solver.getSol();
@@ -530,6 +561,7 @@ int main (int argc, char *argv[]) try
maxSolverSteps,
initialRegularization,
instrumented);
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(xTargetSpace);
solver.solve();
xTargetSpace = solver.getSol();
@@ -577,7 +609,7 @@ int main (int argc, char *argv[]) try
baseTolerance,
instrumented);
- solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(x);
solver.solve();
@@ -616,7 +648,7 @@ int main (int argc, char *argv[]) try
baseTolerance,
instrumented);
- solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(xTargetSpace);
solver.solve();
xTargetSpace = solver.getSol();
@@ -630,6 +662,7 @@ int main (int argc, char *argv[]) try
maxSolverSteps,
initialRegularization,
instrumented);
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(xTargetSpace);
solver.solve();
xTargetSpace = solver.getSol();
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index 8978172ac9fe81b8ce9c8a0b68e5cddc08a9961f..76435d9bf2b9edb6bb05e32908fddc083c2c871d 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -1,4 +1,19 @@
-#define MIXED_SPACE 0
+#ifndef LFE_ORDER
+ #define LFE_ORDER 2
+#endif
+
+#ifndef GFE_ORDER
+ #define GFE_ORDER 2
+#endif
+
+#ifndef MIXED_SPACE
+ #if LFE_ORDER != GFE_ORDER
+ #define MIXED_SPACE 1
+ #else
+ #define MIXED_SPACE 0
+ #endif
+#endif
+
#include <iostream>
#include <fstream>
@@ -63,16 +78,12 @@
#include <dune/solvers/norms/energynorm.hh>
-// grid dimension
-#ifndef WORLD_DIM
-# define WORLD_DIM 3
-#endif
const int dim = WORLD_DIM;
const int targetDim = WORLD_DIM;
-const int displacementOrder = 2;
-const int rotationOrder = 2;
+const int displacementOrder = GFE_ORDER;
+const int rotationOrder = LFE_ORDER;
const int stressFreeDataOrder = 2;
@@ -100,6 +111,9 @@ int main (int argc, char *argv[]) try
#endif
}
+ if (argc < 3)
+ DUNE_THROW(Exception, "Usage: ./film-on-substrate <python path> <parameter file>");
+
// Start Python interpreter
Python::start();
Python::Reference main = Python::import("__main__");
@@ -108,14 +122,13 @@ int main (int argc, char *argv[]) try
//feenableexcept(FE_INVALID);
Python::runStream()
<< std::endl << "import sys"
- << std::endl << "import os"
- << std::endl << "sys.path.append(os.getcwd() + '/../../problems/')"
+ << std::endl << "sys.path.append('" << argv[1] << "')"
<< std::endl;
// parse data file
ParameterTree parameterSet;
- ParameterTreeParser::readINITree(argv[1], parameterSet);
+ ParameterTreeParser::readINITree(argv[2], parameterSet);
ParameterTreeParser::readOptions(argc, argv, parameterSet);
@@ -584,7 +597,7 @@ int main (int argc, char *argv[]) try
baseTolerance,
instrumented);
- solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(x);
solver.solve();
x = solver.getSol();
@@ -604,7 +617,7 @@ int main (int argc, char *argv[]) try
baseTolerance,
instrumented);
- solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(xRBM);
solver.solve();
xRBM = solver.getSol();
@@ -627,6 +640,7 @@ int main (int argc, char *argv[]) try
maxSolverSteps,
initialRegularization,
instrumented);
+ solver.setScaling(parameterSet.get<FieldVector<double,6> >("solverScaling"));
solver.setInitialIterate(xRBM);
solver.solve();
xRBM = solver.getSol();
diff --git a/test/adolctest-scalar-and-vector-mode.cc b/test/adolctest-scalar-and-vector-mode.cc
index 8990660ec830ff719bc9e135e6214f51ff7607dc..036ce5de75b07503afa64f379ec3a535dd5badce 100644
--- a/test/adolctest-scalar-and-vector-mode.cc
+++ b/test/adolctest-scalar-and-vector-mode.cc
@@ -108,6 +108,9 @@ int main (int argc, char *argv[])
parameters["L_c"] = "0.01";
parameters["q"] = "2";
parameters["kappa"] = "1";
+ parameters["b1"] = "1";
+ parameters["b2"] = "1";
+ parameters["b3"] = "1";
//Mixed space
CosseratEnergyLocalStiffness<decltype(compositeBasis), dim,adouble> cosseratEnergyMixed(parameters,
diff --git a/test/adolctest.cc b/test/adolctest.cc
index 2ae115f3157b50917d5e14b27f8c85676da3f665..9a408d0441756b551aef2b1fd5c27b8e47ce9912 100644
--- a/test/adolctest.cc
+++ b/test/adolctest.cc
@@ -490,6 +490,9 @@ int main (int argc, char *argv[]) try
materialParameters["L_c"] = "1";
materialParameters["q"] = "2";
materialParameters["kappa"] = "1";
+ materialParameters["b1"] = "1";
+ materialParameters["b2"] = "1";
+ materialParameters["b3"] = "1";
///////////////////////////////////////////////////////////////////////
// Assemblers for the Euclidean derivatives in an embedding space
diff --git a/test/cosseratenergytest.cc b/test/cosseratenergytest.cc
index 9f4b5bc40e667f1355bf4531f6c25b5d889ecd44..1d41eff2f9fddc017348dc3b04261745ab9faa4f 100644
--- a/test/cosseratenergytest.cc
+++ b/test/cosseratenergytest.cc
@@ -61,6 +61,9 @@ void testEnergy(const GridType* grid, const std::vector<TargetSpace>& coefficien
materialParameters["L_c"] = "0.1";
materialParameters["q"] = "2.5";
materialParameters["kappa"] = "0.1";
+ materialParameters["b1"] = "1";
+ materialParameters["b2"] = "1";
+ materialParameters["b3"] = "1";
typedef Dune::Functions::LagrangeBasis<typename GridType::LeafGridView,1> FEBasis;
FEBasis feBasis(grid->leafGridView());
diff --git a/test/geodesicfeassemblerwrappertest.cc b/test/geodesicfeassemblerwrappertest.cc
index 8030dd8365598d5ff9d91b2686eaa7c10c749612..173a5d7c72df493b6b97e882c19cf035145a4d39 100644
--- a/test/geodesicfeassemblerwrappertest.cc
+++ b/test/geodesicfeassemblerwrappertest.cc
@@ -129,6 +129,9 @@ int main (int argc, char *argv[])
parameters["L_c"] = "0.01";
parameters["q"] = "2";
parameters["kappa"] = "1";
+ parameters["b1"] = "1";
+ parameters["b2"] = "1";
+ parameters["b3"] = "1";
FieldVector<double,dim> values_ = {3e4,2e4,1e4};