From bed86e1a15dbf59a8b33e35f555eac46c9bfe0bf Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Thu, 10 Dec 2020 17:34:04 +0100
Subject: [PATCH] Complete the implementation of the static relaxed
micromorphic model
---
problems/relaxed-micromorphic-torsion.parset | 77 +++
src/relaxed-micromorphic-continuum.cc | 674 ++++++++-----------
2 files changed, 343 insertions(+), 408 deletions(-)
create mode 100644 problems/relaxed-micromorphic-torsion.parset
diff --git a/problems/relaxed-micromorphic-torsion.parset b/problems/relaxed-micromorphic-torsion.parset
new file mode 100644
index 0000000..1147cf7
--- /dev/null
+++ b/problems/relaxed-micromorphic-torsion.parset
@@ -0,0 +1,77 @@
+#############################################
+# Grid parameters
+#############################################
+
+structuredGrid = simplex
+lower = 0 -0.5 -0.5
+upper = 5 0.5 0.5
+elements = 10 2 2
+
+# Number of grid levels
+numLevels = 2
+
+#############################################
+# Solver parameters
+#############################################
+
+timeIntegration = static
+
+# Currently, only 'direct' is supported.
+solverType = direct
+
+# Tolerance of the trust region solver
+#tolerance = 1e-6
+
+# Number of multigrid iterations per trust-region step
+#numIt = 1000
+
+# Number of presmoothing steps
+#nu1 = 3
+
+# Number of postsmoothing steps
+#nu2 = 3
+
+# Number of coarse grid corrections
+#mu = 1
+
+# Number of base solver iterations
+#baseIt = 100
+
+# Tolerance of the multigrid solver
+#mgTolerance = 1e-7
+
+# Tolerance of the base grid solver
+#baseTolerance = 1e-8
+
+############################
+# Material parameters
+############################
+
+energy = stvenantkirchhoff
+
+## For the Wriggers L-shape example
+[materialParameters]
+
+lambda = 1e6
+lambda_h = 1e6
+mu_c = 0
+mu_e = 1e6
+mu_h = 1e6
+
+alpha = 1e3 1e3 1e3
+
+[]
+
+#############################################
+# Boundary values
+#############################################
+
+### Python predicate specifying all Dirichlet grid vertices
+# x is the vertex coordinate
+dirichletVerticesPredicate = "x[0] < 0.0001 or x[0] > 4.9999"
+
+### Neumann values, if needed
+#neumannValues = 0 5e4 0
+
+# Initial iterate (can be a function of x)
+initialIterate = "[[x[0]*1.2, 0, 0], [[0,0,0], [0,0,0], [0,0,0]]]"
diff --git a/src/relaxed-micromorphic-continuum.cc b/src/relaxed-micromorphic-continuum.cc
index c4d10c6..dbca65e 100644
--- a/src/relaxed-micromorphic-continuum.cc
+++ b/src/relaxed-micromorphic-continuum.cc
@@ -3,6 +3,7 @@
#include <dune/common/bitsetvector.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
+#include <dune/common/transpose.hh>
#include <dune/geometry/quadraturerules.hh>
@@ -52,414 +53,258 @@
using namespace Dune;
-
-class RelaxedMicromorphicContinuumAssembler
+/** \brief Return the trace of a matrix
+ if we know that only one row of A has nonzeroes */
+template <class T, int n>
+static T trace(const FieldVector<T,n>& Arow, int row)
{
-#if 0
- // The Lame constants
- double mu_;
- double lambda_;
-
- // The penalty parameter of the DG discretization
- double eta_;
-#endif
- // Neumann boundary force density
- FieldVector<double,2> neumannValue_;
-
-public:
-
- RelaxedMicromorphicContinuumAssembler(const ParameterTree& parameters,
- const FieldVector<double,2> neumannValue)
- : neumannValue_(neumannValue)
- {
-#if 0
- mu_ = parameters.get<double>("mu");
- lambda_ = parameters.get<double>("lambda");
- eta_ = parameters.get<double>("eta");
-#endif
- }
-
- // Compute the stiffness matrix for a single element
- // This local assembler assumes a nonconforming vector-valued basis,
- // and assembles a DG-discretized problem.
- template <class LocalView, class Matrix>
- void assembleElementMatrix(const LocalView& localView, Matrix& elementMatrix) const
- {
- using Element = typename LocalView::Element;
- auto element = localView.element();
- constexpr int dim = Element::dimension;
- auto geometry = element.geometry();
-#if 0
- using Range = typename LocalView::Tree::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
- using Gradient = typename LocalView::Tree::FiniteElement::Traits::LocalBasisType::Traits::JacobianType;
-
- // Get set of shape functions for this element
- const auto& localFiniteElement = localView.tree().finiteElement();
-
- // Set all matrix entries to zero
- elementMatrix.setSize(localFiniteElement.size(),localFiniteElement.size());
- elementMatrix = 0; // Fill the entire matrix with zeros
-
- // Get a quadrature rule
- int order = 2*(dim*localFiniteElement.localBasis().order()-1);
- const auto& quadRule = QuadratureRules<double, dim>::rule(element.type(), order);
+ return Arow[row];
+}
- // Loop over all quadrature points
- for (const auto& qp : quadRule)
+/** \brief Compute the symmetric part of a matrix A, i.e. \f$ \frac 12 (A + A^T) \f$
+ if we know that only one row of A has nonzeroes */
+template <class T, int n>
+static FieldMatrix<T,n,n> sym(const FieldVector<T,n>& Arow, int row)
+{
+ FieldMatrix<T,n,n> result;
+ for (int i=0; i<n; i++)
+ for (int j=0; j<n; j++)
{
- // Position of the current quadrature point in the reference element
- const auto quadPos = qp.position();
-
- // The transposed inverse Jacobian of the map from the reference element to the element
- const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPos);
-
- ///////////////////////////////////////////////////////////////////////////
- // Shape functions - [assume H(div) elements]
- ///////////////////////////////////////////////////////////////////////////
-
- std::vector<Range> values(localFiniteElement.size());
- localFiniteElement.localBasis().evaluateFunction(quadPos, values);
-
- // Gradients of the shape functions on the reference element
- std::vector<Gradient> referenceGradients(localFiniteElement.size());
- localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
-
- // Compute the shape function gradients on the grid element
- std::vector<Gradient> gradients(localFiniteElement.size());
- for (size_t i=0; i<gradients.size(); i++)
- gradients[i] = transformToElementGradient<Gradient,Gradient>(referenceGradients[i], jacobianInverseTransposed);
-
- // geometric weight
- auto factor = qp.weight() * geometry.integrationElement(quadPos);
-
- ///////////////////////////////////
- // Compute the matrix entries
- ///////////////////////////////////
-
- for (size_t i=0; i<localView.size(); i++)
- {
- for (size_t j=0; j<localView.size(); j++)
- {
- // Elastic strains
- FieldMatrix<double,dim,dim> strainU, strainV;
- for (int ii=0; ii<dim; ii++)
- for (int jj=0; jj<dim; jj++)
- {
- strainU[ii][jj] = 0.5 * (gradients[i][ii][jj] + gradients[i][jj][ii]);
- strainV[ii][jj] = 0.5 * (gradients[j][ii][jj] + gradients[j][jj][ii]);
- }
-
- double traceU = 0;
- double traceV = 0;
- for (int ii=0; ii<dim; ii++)
- {
- traceU += strainU[ii][ii];
- traceV += strainV[ii][ii];
- }
-
- double energyDensity = 0;
- for (int ii=0; ii<dim; ii++)
- for (int jj=0; jj<dim; jj++)
- energyDensity += 2 * mu_ * strainU[ii][jj] * strainV[ii][jj];
-
- energyDensity += lambda_ * traceU * traceV;
-
- auto matrixI = localView.tree().localIndex(i);
- auto matrixJ = localView.tree().localIndex(j);
- elementMatrix[matrixI][matrixJ] += energyDensity*factor;
- }
- }
+ T Aij = (i==row) ? Arow[j] : 0;
+ T Aji = (j==row) ? Arow[i] : 0;
+ result[i][j] = 0.5 * (Aij + Aji);
}
-#endif
- }
-
- template <class Intersection, class LocalView, class Matrix>
- void assembleIntersectionMatrix(const Intersection& intersection,
- const LocalView& localViewInside,
- const LocalView& localViewOutside,
- Matrix& matInIn, Matrix& matInOut,
- Matrix& matOutIn, Matrix& matOutOut) const
- {
- // Do nothing -- we are not using a DG discretization yet.
- }
-
- template <class Intersection, class LocalView, class Vector>
- void assembleSurfaceLoad(const Intersection& intersection,
- const LocalView& localView,
- Vector& elementRhs) const
- {
- using Element = typename LocalView::Element;
- constexpr int dim = Element::dimension;
-
- using Range = typename LocalView::Tree::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
-
- auto geometry = intersection.geometry();
- auto geometryInInside = intersection.geometryInInside();
-
- // Get set of shape functions for this element
- const auto& localFiniteElement = localView.tree().finiteElement();
-
- // Set all vector entries to zero
- elementRhs.resize(localFiniteElement.size());
- elementRhs = 0;
-
- // Get a quadrature rule
- int order = 2*(dim*localFiniteElement.localBasis().order()-1);
- const auto& quadRule = QuadratureRules<double, Intersection::mydimension>::rule(intersection.type(), order);
+ return result;
+}
- // Loop over all quadrature points
- for (const auto& qp : quadRule)
+/** \brief Compute the antisymmetric part of a matrix A, i.e. \f$ \frac 12 (A - A^T) \f$
+ if we know that only one row of A has nonzeroes */
+template <class T, int n>
+static FieldMatrix<T,n,n> skew(const FieldVector<T,n>& Arow, int row)
+{
+ FieldMatrix<T,n,n> result;
+ for (int i=0; i<n; i++)
+ for (int j=0; j<n; j++)
{
- // Position of the current quadrature point in the inside element
- const auto quadPosIn = geometryInInside.global(qp.position());
-
- ///////////////////////////////////////////////////////////////////////////
- // Shape functions - [assume H(div) elements]
- ///////////////////////////////////////////////////////////////////////////
-
- std::vector<Range> valuesIn(localFiniteElement.size());
- localFiniteElement.localBasis().evaluateFunction(quadPosIn, valuesIn);
-
- // geometric weight
- auto factor = qp.weight() * geometry.integrationElement(qp.position());
-
- ///////////////////////////////////
- // Compute the matrix entries
- ///////////////////////////////////
- for (size_t i=0; i<localView.size(); i++)
- {
- auto matrixI = localView.tree().localIndex(i);
- elementRhs[matrixI] += valuesIn[i] * neumannValue_ * factor;
- }
-
+ T Aij = (i==row) ? Arow[j] : 0;
+ T Aji = (j==row) ? Arow[i] : 0;
+ result[i][j] = 0.5 * (Aij - Aji);
}
- }
-
-};
+ return result;
+}
-// Get the occupation pattern of the stiffness matrix
-template <class Basis, class MatrixType>
-void setOccupationPattern(const Basis& basis, MatrixType& matrix)
+/** \brief Compute the deviator of a matrix A */
+template <class T, int n>
+static FieldMatrix<T,n,n> dev(const FieldMatrix<T,n,n>& A)
{
- // MatrixIndexSets store the occupation pattern of a sparse matrix.
- // They are not particularly efficient, but simple to use.
- std::array<std::array<MatrixIndexSet, 2>, 2> nb;
-
- // Set sizes of the 2x2 submatrices
- for (size_t i=0; i<2; i++)
- for (size_t j=0; j<2; j++)
- nb[i][j].resize(basis.size({i}), basis.size({j}));
-
- // A view on the FE basis on a single element
- auto localView = basis.localView();
-
- // Loop over all leaf elements
- for(const auto& element : elements(basis.gridView()))
- {
- // Bind the local view to the current element
- localView.bind(element);
-
- // Add element stiffness matrix onto the global stiffness matrix
- for (size_t i=0; i<localView.size(); i++) {
-
- // Global index of the i-th local degree of freedom of the current element
- auto row = localView.index(i);
-
- for (size_t j=0; j<localView.size(); j++ ) {
+ FieldMatrix<T,n,n> result = A;
+ T trace(0);
+ for (int i=0; i<n; i++)
+ trace += A[i][i];
+
+ for (int i=0; i<n; i++)
+ result[i][i] -= trace / n;
+ return result;
+}
- // Global index of the j-th local degree of freedom of the current element
- auto col = localView.index(j);
+/** \brief The Frobenius (i.e., componentwise) product of two matrices */
+template <class T, int n>
+static T frobeniusProduct(const FieldMatrix<T,n,n>& A, const FieldMatrix<T,n,n>& B)
+{
+ T result(0.0);
- nb[row[0]][col[0]].add(row[1],col[1]);
- }
- }
- }
+ for (int i=0; i<n; i++)
+ for (int j=0; j<n; j++)
+ result += A[i][j] * B[i][j];
- // Give the matrix the occupation pattern we want.
- using namespace Indices;
- nb[0][0].exportIdx(matrix[_0][_0]);
- nb[0][1].exportIdx(matrix[_0][_1]);
- nb[1][0].exportIdx(matrix[_1][_0]);
- nb[1][1].exportIdx(matrix[_1][_1]);
+ return result;
}
-template<class Matrix, class MultiIndex>
-decltype(auto) matrixEntry(
- Matrix& matrix, const MultiIndex& row, const MultiIndex& col)
+// \begin{equation*}
+// (\Curl P)_i
+// \colonequals
+// \Big( \frac{\partial P_3}{\partial x_2} - \frac{\partial P_2}{\partial x_3}, \quad
+// \frac{\partial P_1}{\partial x_3} - \frac{\partial P_3}{\partial x_1}, \quad
+// \frac{\partial P_2}{\partial x_1} - \frac{\partial P_1}{\partial x_2} \Big).
+// \end{equation*}
+template <class T, int n>
+static FieldVector<T,n> curl(const FieldMatrix<T,n,n>& dP)
{
- using namespace Indices;
- if ((row[0]==0) and (col[0]==0))
- return matrix[_0][_0][row[1]][col[1]][row[2]][col[2]];
- if ((row[0]==0) and (col[0]==1))
- return matrix[_0][_1][row[1]][col[1]][row[2]][0];
- if ((row[0]==1) and (col[0]==0))
- return matrix[_1][_0][row[1]][col[1]][0][col[2]];
- return matrix[_1][_1][row[1]][col[1]][0][0];
+ return {dP[2][1] - dP[1][2], dP[0][2] - dP[2][0], dP[1][0] - dP[0][1]};
}
-
-/** \brief Assemble the stiffness matrix on the given grid view */
-template <class Basis, class LocalAssembler, class Matrix, class Vector>
-void assembleStiffnessMatrix(const Basis& basis,
- const LocalAssembler& localAssembler,
- const BoundaryPatch<typename Basis::GridView>& neumannBoundary,
- Matrix& matrix,
- Vector& rhs)
+/** \brief Local mass assembler for dune-functions basis
+ *
+ * We assume the mass matrix to be symmetric and hence ansatz and trial space to be equal!
+ * The implementation allows an arbitrary dune-functions compatible basis, as long as there
+ * is an ISTLBackend available for the resulting matrix.
+ */
+class RelaxedMicromorphicContinuumAssembler
{
- auto gridView = basis.gridView();
+ double mu_e_;
+ double mu_c_;
+ double lambda_;
- setOccupationPattern(basis, matrix);
+ double mu_h_;
+ double lambda_h_;
+ std::array<double,3> alpha_;
+public:
- // Unique element indices for visiting every intersection only once
-// using ElementMapper = MultipleCodimMultipleGeomTypeMapper<typename Basis::GridView>;
-// ElementMapper elementMapper(gridView, mcmgElementLayout());
+ RelaxedMicromorphicContinuumAssembler(const ParameterTree& parameters)
+ {
+ mu_e_ = parameters.get<double>("mu_e");
+ mu_c_ = parameters.get<double>("mu_c");
+ lambda_ = parameters.get<double>("lambda");
- // Set all entries to zero
- matrix = 0;
+ mu_h_ = parameters.get<double>("mu_h");
+ lambda_h_ = parameters.get<double>("lambda_h");
+ alpha_ = parameters.get<std::array<double,3> >("alpha");
+ }
- // Set rhs to correct length
- auto rhsBackend = Functions::istlVectorBackend(rhs);
- rhsBackend.resize(basis);
+ template<class Element, class LocalMatrix, class LocalView>
+ void operator()(const Element& element, LocalMatrix& localMatrix, const LocalView& trialLocalView, const LocalView& ansatzLocalView) const
+ {
+ using namespace Dune::Indices;
- // Set all entries to zero
- rhs = 0;
+ // the mass matrix operator assumes trial == ansatz
+ if (&trialLocalView != &ansatzLocalView)
+ DUNE_THROW(Dune::Exception,"The RelaxedMicromorphicContinuumAssembler assumes equal ansatz and trial bases");
- // A loop over all elements of the grid
- auto localView = basis.localView();
+ // matrix was already resized before
+ localMatrix = 0.;
- for (const auto& element : elements(gridView))
- {
- localView.bind(element);
+ const auto& geometry = element.geometry();
+ constexpr int gridDim = LocalView::GridView::dimension;
- // Now let's get the element stiffness matrix
- // A dense matrix is used for the element stiffness matrix
- Dune::Matrix<double> elementMatrix;
- localAssembler.assembleElementMatrix(localView, elementMatrix);
+ const auto& displacementFiniteElement = trialLocalView.tree().child(_0,0).finiteElement();
+ const auto& localDisplacementBasis = displacementFiniteElement.localBasis();
+ using DisplacementJacobianType = typename std::decay_t<decltype(localDisplacementBasis)>::Traits::JacobianType;
- // Add element stiffness matrix onto the global stiffness matrix
- for (size_t i=0; i<elementMatrix.N(); i++)
- {
- // The global index of the i-th local degree of freedom of the element
- auto row = localView.index(i);
+ const auto& distortionFiniteElement = trialLocalView.tree().child(_1,0).finiteElement();
+ const auto& localDistortionBasis = distortionFiniteElement.localBasis();
+ using DistortionRangeType = typename std::decay_t<decltype(localDistortionBasis)>::Traits::RangeType;
+ using DistortionJacobianType = typename std::decay_t<decltype(localDistortionBasis)>::Traits::JacobianType;
- for (size_t j=0; j<elementMatrix.M(); j++ )
- {
- // The global index of the j-th local degree of freedom of the element
- auto col = localView.index(j);
-#if 0 // If matrix is a MultiTypeBlockMatrix
- matrixEntry(matrix, row, col) += elementMatrix[i][j];
-#else
- matrix[row][col] += elementMatrix[i][j];
-#endif
- }
- }
- }
+ std::vector<DisplacementJacobianType> displacementJacobian(localDisplacementBasis.size());
-#if 0
- auto localViewOutside = basis.localView();
+ std::vector<DistortionRangeType> microDistortion(localDistortionBasis.size());
+ std::vector<DistortionJacobianType> microDistortionJacobian(localDistortionBasis.size());
+ std::vector<FieldVector<double,gridDim> > microDistortionCurl(localDistortionBasis.size());
- // A loop over all intersections of the grid
- for (const auto& inside : elements(gridView))
- {
- localView.bind(inside);
+ const int quadOrder = gridDim; // TODO: More appropriate order?
+ const auto& quad = QuadratureRules<double, gridDim>::rule(element.type(), quadOrder);
- for (const auto& intersection : intersections(gridView, inside))
+ for( const auto& qp : quad )
{
- if (intersection.neighbor())
- {
- // Visit every intersection only once
- if (elementMapper.index(inside) > elementMapper.index(intersection.outside()))
- continue;
+ const auto quadPos = qp.position();
+ const auto integrationElement = geometry.integrationElement(quadPos);
+ const auto geometryJacobianIT = geometry.jacobianInverseTransposed(quadPos);
- localViewOutside.bind(intersection.outside());
+ localDisplacementBasis.evaluateJacobian(quadPos, displacementJacobian);
- Dune::Matrix<double> matInIn;
- Dune::Matrix<double> matInOut;
- Dune::Matrix<double> matOutIn;
- Dune::Matrix<double> matOutOut;
+ for (std::size_t i=0; i<displacementJacobian.size(); i++)
+ displacementJacobian[i] = displacementJacobian[i] * transpose(geometryJacobianIT);
- localAssembler.assembleIntersectionMatrix(intersection,
- localView, localViewOutside,
- matInIn, matInOut, matOutIn, matOutOut);
+ localDistortionBasis.evaluateFunction(quadPos, microDistortion);
- for(size_t i=0; i<matInIn.N(); i++)
- {
- // The global index of the i-th degree of freedom of the element
- auto row = localView.index(i);
+ localDistortionBasis.evaluateJacobian(quadPos, microDistortionJacobian);
- for (size_t j=0; j<matInIn.M(); j++ )
- {
- // The global index of the j-th degree of freedom of the element
- auto col = localView.index(j);
- matrix[row][col] += matInIn[i][j];
- }
- }
+ for (std::size_t i=0; i<microDistortionJacobian.size(); i++)
+ {
+ microDistortionJacobian[i] = microDistortionJacobian[i] * transpose(geometryJacobianIT);
+ microDistortionCurl[i] = curl(microDistortionJacobian[i]);
+ }
- for(size_t i=0; i<matInOut.N(); i++)
+ // The displacement coupling with itself
+ for (std::size_t i=0; i<localDisplacementBasis.size(); i++)
+ {
+ for (std::size_t iDir=0; iDir<gridDim; iDir++)
{
- // The global index of the i-th degree of freedom of the element
- auto row = localView.index(i);
+ auto rowIndex = trialLocalView.tree().child(_0,iDir).localIndex(i);
- for (size_t j=0; j<matInOut.M(); j++ )
+ for (std::size_t j=0; j<localDisplacementBasis.size(); j++)
{
- // The global index of the j-th degree of freedom of the element
- auto col = localViewOutside.index(j);
- matrix[row][col] += matInOut[i][j];
+ for (std::size_t jDir=0; jDir<gridDim; jDir++)
+ {
+ auto colIndex = trialLocalView.tree().child(_0,jDir).localIndex(j);
+
+ auto value = mu_e_ * frobeniusProduct(sym(displacementJacobian[i][0], iDir),
+ sym(displacementJacobian[j][0], jDir))
+ + mu_c_ * frobeniusProduct(skew(displacementJacobian[i][0], iDir),
+ skew(displacementJacobian[j][0], jDir))
+ + 0.5 * lambda_ * trace(displacementJacobian[i][0], iDir) * trace(displacementJacobian[j][0], jDir);
+ localMatrix[rowIndex][colIndex] += value * qp.weight() * integrationElement;
+ }
}
}
+ }
- for(size_t i=0; i<matOutIn.N(); i++)
+ // The displacement coupling with the micro-distortion
+ for (std::size_t i=0; i<localDisplacementBasis.size(); i++)
+ {
+ for (std::size_t iDir=0; iDir<gridDim; iDir++)
{
- // The global index of the i-th degree of freedom of the element
- auto row = localViewOutside.index(i);
+ auto rowIndex = trialLocalView.tree().child(_0,iDir).localIndex(i);
- for (size_t j=0; j<matOutIn.M(); j++ )
+ for (std::size_t j=0; j<localDistortionBasis.size(); j++)
{
- // The global index of the j-th degree of freedom of the element
- auto col = localView.index(j);
- matrix[row][col] += matOutIn[i][j];
- }
- }
+ for (std::size_t jDir=0; jDir<gridDim; jDir++)
+ {
+ auto colIndex = trialLocalView.tree().child(_1,jDir).localIndex(j);
- for (size_t i=0; i<matOutOut.N(); i++)
- {
- // The global index of the i-th degree of freedom of the element
- auto row = localViewOutside.index(i);
+ auto value = -mu_e_ * frobeniusProduct(sym(displacementJacobian[i][0], iDir),
+ sym(microDistortion[j], jDir))
+ - mu_c_ * frobeniusProduct(skew(displacementJacobian[i][0], iDir),
+ skew(microDistortion[j], jDir))
+ - 0.5 * lambda_ * trace(displacementJacobian[i][0], iDir) * trace(microDistortion[j], jDir);
- for (size_t j=0; j<matOutOut.M(); j++ )
- {
- // The global index of the j-th degree of freedom of the element
- auto col = localViewOutside.index(j);
- matrix[row][col] += matOutOut[i][j];
+ localMatrix[rowIndex][colIndex] += value * qp.weight() * integrationElement;
+ localMatrix[colIndex][rowIndex] += value * qp.weight() * integrationElement;
+ }
}
}
+
}
- else if (neumannBoundary.contains(intersection))
+
+ // The micro-distortion coupling with itself
+ for (std::size_t i=0; i<localDistortionBasis.size(); i++)
{
- Vector elementRhs(localView.size());
- elementRhs = 0;
+ for (std::size_t iDir=0; iDir<gridDim; iDir++)
+ {
+ auto rowIndex = trialLocalView.tree().child(_1,iDir).localIndex(i);
- localAssembler.assembleSurfaceLoad(intersection,
- localView,
- elementRhs);
+ for (std::size_t j=0; j<localDistortionBasis.size(); j++)
+ {
+ for (std::size_t jDir=0; jDir<gridDim; jDir++)
+ {
+ auto colIndex = trialLocalView.tree().child(_1,jDir).localIndex(j);
+ auto value = mu_e_ * frobeniusProduct(sym(microDistortion[i], iDir), sym(microDistortion[j], jDir))
+ + mu_c_ * frobeniusProduct(skew(microDistortion[i], iDir), skew(microDistortion[j], jDir))
+ + 0.5 * lambda_ * trace(microDistortion[i], iDir) * trace(microDistortion[j], jDir);
- for (size_t i=0; i<localView.size(); i++)
- {
- // The global index of the i-th degree of freedom of the element
- auto row = localView.index(i);
+ value += mu_h_ * frobeniusProduct(sym(microDistortion[i], iDir), sym(microDistortion[j], jDir))
+ + 0.5 * lambda_h_ * trace(microDistortion[i], iDir) * trace(microDistortion[j], jDir);
+
+ value += 0.5 * alpha_[0] * frobeniusProduct(dev(sym(microDistortionCurl[i], iDir)),
+ dev(sym(microDistortionCurl[j], jDir)))
+ + 0.5 * alpha_[1] * frobeniusProduct(skew(microDistortionCurl[i], iDir),
+ skew(microDistortionCurl[j], jDir))
+ + 0.5 * alpha_[2] * trace(microDistortionCurl[i], iDir) * trace(microDistortionCurl[j], jDir);
- rhs[row] += elementRhs[i];
+ localMatrix[rowIndex][colIndex] += value * qp.weight() * integrationElement;
+ }
+ }
}
}
}
}
-#endif
-}
+
+};
// The grid dimension
-const int dim = 2;
+const int dim = 3;
int main (int argc, char *argv[]) try
{
@@ -491,13 +336,6 @@ int main (int argc, char *argv[]) try
// read solver settings
const std::string solverType = parameterSet["solverType"];
- // Multigrid settings, if the solver is a multigrid one
- const int maxIterations = parameterSet.get<int>("maxIterations");
- const int nu1 = parameterSet.get<int>("nu1");
- const int nu2 = parameterSet.get<int>("nu2");
- const int mu = parameterSet.get<int>("mu");
- const double tolerance = parameterSet.get<double>("tolerance");
-
///////////////////////////////
// Generate the grid
///////////////////////////////
@@ -548,7 +386,6 @@ int main (int argc, char *argv[]) try
using namespace Indices;
auto displacementBasis = Functions::subspaceBasis(basis, _0);
- auto microDistortionBasis = Functions::subspaceBasis(basis, _1);
///////////////////////////////////////////
@@ -571,50 +408,12 @@ int main (int argc, char *argv[]) try
BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
- ///////////////////////////////////////////
- // Construct Neumann boundary
- ///////////////////////////////////////////
- BitSetVector<1> neumannVertices(gridView.size(dim), false);
-
- // Make Python function that computes which vertices are on the Dirichlet boundary,
- // based on the vertex positions.
- lambda = std::string("lambda x: (") + parameterSet.get<std::string>("neumannVerticesPredicate") + std::string(")");
- auto pythonNeumannVertices = Python::make_function<bool>(Python::evaluate(lambda));
-
- for (auto&& vertex : vertices(gridView))
- {
- bool isNeumann = pythonNeumannVertices(vertex.geometry().corner(0));
- neumannVertices[indexSet.index(vertex)] = isNeumann;
- }
-
- BoundaryPatch<GridView> neumannBoundary(gridView, neumannVertices);
-
- BitSetVector<1> dirichletNodes(basis.size(), false);
- constructBoundaryDofs(dirichletBoundary,basis,dirichletNodes);
+ BitSetVector<1> dirichletDofs(basis.size(), false);
+ constructBoundaryDofs(dirichletBoundary,basis,dirichletDofs);
// Vector and matrix types
-#if 0
- // Right now there is no real need for MultiTypeBlockVector, because both displacement
- // and micro-distortion use the same vector type. However, I suspect we will have to
- // Raviart-Thomas elements for the displacement eventually, so let's prepare for the
- // future and use MultiTypeBlockVector right away.
- using DisplacementVector = BlockVector<FieldVector<double,dim> >;
- using microDistortionVector = BlockVector<FieldVector<double,dim> >;
-
- using Vector = MultiTypeBlockVector<DisplacementVector, microDistortionVector>;
-
- using Matrix00 = BCRSMatrix<FieldMatrix<double,dim,dim>>;
- using Matrix01 = BCRSMatrix<FieldMatrix<double,dim,dim>>;
- using Matrix10 = BCRSMatrix<FieldMatrix<double,dim,dim>>;
- using Matrix11 = BCRSMatrix<FieldMatrix<double,dim,dim>>;
- using MatrixRow0 = MultiTypeBlockVector<Matrix00, Matrix01>;
- using MatrixRow1 = MultiTypeBlockVector<Matrix10, Matrix11>;
- using Matrix = MultiTypeBlockMatrix<MatrixRow0,MatrixRow1>;
-#else
using Vector = BlockVector<double>;
using Matrix = BCRSMatrix<double>;
-#endif
-
////////////////////////////
// Initial iterate
@@ -633,22 +432,27 @@ int main (int argc, char *argv[]) try
rhs = 0;
// Assemble elasticity problem
- RelaxedMicromorphicContinuumAssembler localAssembler(parameterSet.sub("materialParameters"),
- parameterSet.get<FieldVector<double,dim> >("neumannValue"));
+ RelaxedMicromorphicContinuumAssembler localAssembler(parameterSet.sub("materialParameters"));
Matrix stiffnessMatrix;
- assembleStiffnessMatrix(basis, localAssembler, neumannBoundary, stiffnessMatrix, rhs);
+ auto assembler = Fufem::duneFunctionsOperatorAssembler(basis, basis);
+
+ Timer watch;
+ assembler.assembleBulk(Fufem::istlMatrixBackend(stiffnessMatrix), localAssembler);
+ std::cout << "Assembly took " << watch.elapsed() << " seconds" << std::endl;
-#if 0
///////////////////////////////////////////
// Modify Dirichlet rows
///////////////////////////////////////////
for (size_t i=0; i<stiffnessMatrix.N(); i++)
- if (dirichletNodes[i][0])
+ if (dirichletDofs[i][0])
+ {
for (auto&& [entry, j] : sparseRange(stiffnessMatrix[i]))
entry = (i==j) ? 1.0 : 0.0;
-#endif
+
+ rhs[i] = x[i];
+ }
/////////////////////////////////
// Assemble the mass matrix
@@ -664,15 +468,22 @@ int main (int argc, char *argv[]) try
// Create a solver
/////////////////////////////
+ std::shared_ptr<Solvers::LinearSolver<Matrix,Vector> > solver;
+
if (solverType == "direct")
{
- Solvers::UMFPackSolver<Matrix,Vector> solver(stiffnessMatrix, x, rhs);
-
- solver.solve();
+ solver = std::make_shared<Solvers::UMFPackSolver<Matrix,Vector> >();
}
#if 0
else if (solverType == "standardMG")
{
+ // Multigrid settings, if the solver is a multigrid one
+ const int maxIterations = parameterSet.get<int>("maxIterations");
+ const int nu1 = parameterSet.get<int>("nu1");
+ const int nu2 = parameterSet.get<int>("nu2");
+ const int mu = parameterSet.get<int>("mu");
+ const double tolerance = parameterSet.get<double>("tolerance");
+
// Make pre and postsmoothers
auto presmoother = Solvers::BlockGSStepFactory<Matrix, Vector>::create(
Solvers::BlockGS::LocalSolvers::direct(0.0));
@@ -725,15 +536,62 @@ int main (int argc, char *argv[]) try
else
DUNE_THROW(NotImplemented, "Unknown solver type!");
- ///////////////////////////
- // Write result to file
- ///////////////////////////
+ if (parameterSet["timeIntegration"] == "static")
+ {
+ solver->setProblem(stiffnessMatrix, x, rhs);
+ solver->solve();
+
+ ///////////////////////////
+ // Write result to file
+ ///////////////////////////
+
+ auto displacementFunction = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim> >(displacementBasis, x);
+
+ auto director0Function = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim> >(Functions::subspaceBasis(basis, _1, 0), x);
+ auto director1Function = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim> >(Functions::subspaceBasis(basis, _1, 1), x);
+ auto director2Function = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim> >(Functions::subspaceBasis(basis, _1, 2), x);
+
+ VTKWriter<GridView> vtkWriter(grid->leafGridView(), VTK::nonconforming);
+ vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
+ vtkWriter.addVertexData(director0Function, VTK::FieldInfo("director0", VTK::FieldInfo::Type::vector, dim));
+ vtkWriter.addVertexData(director1Function, VTK::FieldInfo("director1", VTK::FieldInfo::Type::vector, dim));
+ vtkWriter.addVertexData(director2Function, VTK::FieldInfo("director2", VTK::FieldInfo::Type::vector, dim));
+ vtkWriter.write("relaxed-micromorphic-continuum-result");
+ }
+ else if (parameterSet["timeIntegration"] == "implicitMidpoint")
+ {
+ double timeStep = parameterSet.get<double>("timeStep");
+
+ // Set up the matrix for the increment problem
+ auto matrix = massMatrix;
+ massMatrix.axpy(power(timeStep,2) * 0.25, stiffnessMatrix);
+#if 0
+ for (int i=0; i<numTimeSteps; i++)
+ {
+ // First half-step: Solve for the new configuration
+ auto rhs = M * oldConfiguration
+ + timeStep * 0.5 * (oldMomentum - timeStep * 0.5 * stiffnessMatrix * oldConfiguration - timeStep * b);
+
+ solver->setProblem(matrix, newConfiguration, rhs);
+ solver->solve();
- auto displacementFunction = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim> >(displacementBasis, x);
+ // Second half-step: Update momenta
+ newMomentum = momentum - timeStep * 0.5 * stiffnessMatrix * (oldConfiguration + newConfiguration) - 2*b;
+
+ ///////////////////////////
+ // Write result to file
+ ///////////////////////////
+
+ auto displacementFunction = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim> >(displacementBasis, x);
+
+ VTKWriter<GridView> vtkWriter(grid->leafGridView(), VTK::nonconforming);
+ vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
+ vtkWriter.write("relaxed-micromorphic-continuum-result");
+ }
+#endif
+ } else
+ DUNE_THROW(NotImplemented, "Time integration method '" << parameterSet["timeIntegration"] << "' not implemented!");
- VTKWriter<GridView> vtkWriter(grid->leafGridView(), VTK::nonconforming);
- vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
- vtkWriter.write("relaxed-micromorphic-continuum-result");
} catch (Exception& e)
{
std::cout << e.what() << std::endl;
--
GitLab