Skip to content
Snippets Groups Projects
Commit bed86e1a authored by oliver.sander_at_tu-dresden.de's avatar oliver.sander_at_tu-dresden.de
Browse files

Complete the implementation of the static relaxed micromorphic model

parent 06205905
No related branches found
No related tags found
No related merge requests found
Pipeline #33815 failed
#############################################
# 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]]]"
......@@ -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;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment