Skip to content

GitLab

Commit 2fde2da4 authored by Patrick Jaap's avatar Patrick Jaap
Browse files

Merge branch 'cleanup/whitespaces' into 'master' 

Whitespace cleanup

See merge request !29
parents cc26b05c be66c63f
Hide whitespace changes
Inline Side-by-side
dune/elasticity/assemblers/mooneyrivlinfunctionalassembler.hh
View file @ 2fde2da4
......@@ -15,7 +15,7 @@
/** \brief Local assembler for the linearization of a nonlinear Mooney-Rivlin energy at a displacement u
*
* \f[
* W(u)= a tr(E(u)) + b tr(E(u))^2 + c ||E(u)||^2 + d J^2 + e J^-k, k>=2
* W(u)= a tr(E(u)) + b tr(E(u))^2 + c ||E(u)||^2 + d J^2 + e J^-k, k>=2
* \f]
*
* where
......@@ -24,9 +24,9 @@
* - \f$J\f$ the determinant of the deformation gradient
* - \f$\a\f$,...,\f$\e\f$ material parameters
* - \f$k\f$ controlls orientation preservation (k \approx \floor(1/(1-2nu) -1))
*/
*/
template <class GridType, class TrialLocalFE>
class MooneyRivlinFunctionalAssembler :
class MooneyRivlinFunctionalAssembler :
public LocalFunctionalAssembler<GridType,TrialLocalFE, Dune::FieldVector<typename GridType::ctype ,GridType::dimension> >
{
......@@ -58,7 +58,7 @@ public:
localVector = 0.0;
// get quadrature rule (should depend on k?)
const int order = (element.type().isSimplex())
const int order = (element.type().isSimplex())
? 6*(tFE.localBasis().order())
: 6*(tFE.localBasis().order()*dim);
const auto& quad = QuadratureRuleCache<ctype, dim>::rule(element.type(), order, IsRefinedLocalFiniteElement<TrialLocalFE>::value(tFE) );
......@@ -112,7 +112,7 @@ public:
linStrain[i][k] = Dune::Elasticity::linearisedStrain(deformationGradient, localConfGrad);
}
// collect terms
// collect terms
FMatrix fu(0);
// the linearization of the trace is just deformation gradient
......
dune/elasticity/assemblers/mooneyrivlinoperatorassembler.hh
View file @ 2fde2da4
......@@ -18,7 +18,7 @@
/** \brief Assemble the second derivative of a Mooney-Rivlin material.
*
* \f[
* W(u)= a tr(E(u)) + b tr(E(u))^2 + c ||E(u)||^2 + d J^2 + e J^-k, k>=2
* W(u)= a tr(E(u)) + b tr(E(u))^2 + c ||E(u)||^2 + d J^2 + e J^-k, k>=2
* \f]
*
* where
......@@ -27,9 +27,9 @@
* - \f$J\f$ the determinant of the deformation gradient
* - \f$\a\f$,...,\f$\e\f$ material parameters
* - \f$k\f$ controlls orientation preservation (k \approx \floor(1/(1-2nu) -1))
*/
*/
template <class GridType, class TrialLocalFE, class AnsatzLocalFE>
class MooneyRivlinOperatorAssembler
class MooneyRivlinOperatorAssembler
: public LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE, Dune::FieldMatrix<typename GridType::ctype ,GridType::dimension,GridType::dimension> >
{
static constexpr int dim = GridType::dimension;
......@@ -67,7 +67,7 @@ public:
localMatrix = 0.0;
// get quadrature rule (should depend on k?)
const int order = (element.type().isSimplex())
const int order = (element.type().isSimplex())
? 6*(tFE.localBasis().order())
: 6*(tFE.localBasis().order()*dim);
const auto& quad = QuadratureRuleCache<ctype, dim>::rule(element.type(), order, IsRefinedLocalFiniteElement<TrialLocalFE>::value(tFE));
......@@ -131,7 +131,7 @@ public:
for (int row = 0; row < rows; ++row)
for (int col = 0; col<cols; ++col) {
// second derivative of the trace terms
auto hessTrace = z*(a_ + 2*b_*strain.trace()) * ((gradients[row]*gradients[col]));
Dune::MatrixVector::addToDiagonal(localMatrix[row][col], hessTrace);
......@@ -180,14 +180,14 @@ private:
std::shared_ptr<GridFunction> displacement_;
/** \brief Compute the matrix entry that arises from the second derivative of the determinant.
* Note that it is assumed the defGrad is the deformation gradient and not the displacement gradient
/** \brief Compute the matrix entry that arises from the second derivative of the determinant.
* Note that it is assumed the defGrad is the deformation gradient and not the displacement gradient
*/
auto hessDefDet(const Dune::FieldMatrix<ctype,dim,dim>& defGrad,
const Dune::FieldVector<ctype,dim>& testGrad, const Dune::FieldVector<ctype,dim>& ansatzGrad) const
{
Dune::FieldMatrix<ctype, dim, dim> res(0);
if (dim==2) {
res[0][1] = testGrad[0]*ansatzGrad[1]-ansatzGrad[0]*testGrad[1];
res[1][0] = -res[0][1];
......@@ -195,7 +195,7 @@ private:
}
//compute the cross product
Dune::FieldVector<ctype,dim> cross(0);
Dune::FieldVector<ctype,dim> cross(0);
for (int k=0; k<dim; k++)
cross[k]=testGrad[(k+1)%dim]*ansatzGrad[(k+2)%dim] -testGrad[(k+2)%dim]*ansatzGrad[(k+1)%dim];
......@@ -206,7 +206,7 @@ private:
res[i][j] = (cross*defGrad[k])*std::pow(-1,k);
res[j][i] = -res[i][j];
}
return res;
}
};
......
dune/elasticity/assemblers/neohookefunctionalassembler.hh
View file @ 2fde2da4
......@@ -25,10 +25,10 @@
*
* The linearization of this energy at \f$ u\f$ as a functional applied to a test function \f$ v\f$] reads
*
*
*
* Laursen:
* \f[
* l_u(v) := (\frac{\lambda}2 J - \frac{\frac{\lambda}2 + \mu}{J})\frac{\partial J}{\partial u}v +
* l_u(v) := (\frac{\lambda}2 J - \frac{\frac{\lambda}2 + \mu}{J})\frac{\partial J}{\partial u}v +
* \mu\frac{\partial tr(E(u))}{\partial u} v
* \f]
*
......@@ -43,9 +43,9 @@
* - \f$tr \f$: the trace operator
* - \f$J\f$ the determinant of the deformation gradient
* - \f$\lambda\f$,\f$\mu\f$ material parameters (first lame and shear modulus)
*/
*/
template <class GridType, class TrialLocalFE>
class NeoHookeFunctionalAssembler :
class NeoHookeFunctionalAssembler :
public LocalFunctionalAssembler<GridType, TrialLocalFE, Dune::FieldVector<typename GridType::ctype ,GridType::dimension> >
{
......@@ -78,8 +78,8 @@ public:
{
typedef Dune::FieldMatrix<ctype,dim,dim> FMdimdim;
typedef typename Element::Geometry::JacobianInverseTransposed JacInvTransType;
typedef typename TrialLocalFE::Traits::LocalBasisType::Traits::JacobianType JacobianType;
typedef typename GridType::template Codim<0>::Geometry::LocalCoordinate LocalCoordinate;
typedef typename TrialLocalFE::Traits::LocalBasisType::Traits::JacobianType JacobianType;
typedef typename GridType::template Codim<0>::Geometry::LocalCoordinate LocalCoordinate;
int rows = tFE.localBasis().size();
......@@ -138,12 +138,12 @@ public:
// make deformation gradient out of the discplacement
for (int i=0;i<dim;i++)
localConfGrad[i][i] += 1;
localConfGrad[i][i] += 1;
// evaluate the derminante of the deformation gradient
const ctype J = localConfGrad.determinant();
// collect terms
// collect terms
FMdimdim fu(0);
#ifdef LAURSEN
......@@ -159,7 +159,7 @@ public:
ctype z = quad[pt].weight()*integrationElement;
// add vector entries
for(int row=0; row<rows; row++)
for(int row=0; row<rows; row++)
fu.usmv(z,gradients[row],localVector[row]);
}
}
......
dune/elasticity/assemblers/neohookeoperatorassembler.hh
View file @ 2fde2da4
......@@ -50,14 +50,14 @@
* - \f$\lambda\f$,\f$\mu\f$ material parameters (first lame and shear modulus)
*/
template <class GridType, class TrialLocalFE, class AnsatzLocalFE>
class NeoHookeOperatorAssembler
class NeoHookeOperatorAssembler
: public LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE, Dune::FieldMatrix<typename GridType::ctype ,GridType::dimension,GridType::dimension> >
{
static const int dim = GridType::dimension;
typedef typename GridType::ctype ctype;
typedef VirtualGridFunction<GridType, Dune::FieldVector<ctype,GridType::dimension> > GridFunction;
public:
......@@ -95,7 +95,7 @@ public:
localMatrix = 0.0;
// TODO get proper quadrature rule
const int order = (element.type().isSimplex())
const int order = (element.type().isSimplex())
? 4*(tFE.localBasis().order())
: 4*(tFE.localBasis().order()*dim);
const Dune::QuadratureRule<ctype, dim>& quad = QuadratureRuleCache<ctype, dim>::rule(element.type(), order, IsRefinedLocalFiniteElement<TrialLocalFE>::value(tFE) );
......@@ -144,7 +144,7 @@ public:
// compute deformation gradient of the configuration
for (int i=0;i<dim;i++)
localConfGrad[i][i] += 1;
localConfGrad[i][i] += 1;
// evaluate the derminante of the deformation gradient
const ctype J = localConfGrad.determinant();
......@@ -166,18 +166,18 @@ public:
for (int row=0; row<rows; row++)
for (int col=0; col<cols; col++) {
// second derivative of the determinat term
// second derivative of the determinat term
localMatrix[row][col].axpy(z*coeff,hessDefDet(localConfGrad,gradients[row],gradients[col]));
FVdim rowTerm(0),colTerm(0);
linDefDet.mv(gradients[row],rowTerm);
linDefDet.mv(gradients[col],colTerm);
linDefDet.mv(gradients[row],rowTerm);
linDefDet.mv(gradients[col],colTerm);
// derivative of the coefficient term
for (int rcomp=0; rcomp<dim; rcomp++)
localMatrix[row][col][rcomp].axpy(coeff2*z*rowTerm[rcomp],colTerm);
// second derivative of the trace term
StaticMatrix::addToDiagonal(localMatrix[row][col],z*mu_*((gradients[row]*gradients[col])));
......@@ -202,14 +202,14 @@ private:
std::shared_ptr<GridFunction> displacement_;
/** \brief Compute the matrix entry that arises from the second derivative of the determinant.
* Note that it is assumed the defGrad is the deformation gradient and not the displacement gradient
/** \brief Compute the matrix entry that arises from the second derivative of the determinant.
* Note that it is assumed the defGrad is the deformation gradient and not the displacement gradient
*/
Dune::FieldMatrix<ctype,dim,dim> hessDefDet(const Dune::FieldMatrix<ctype,dim,dim>& defGrad,
const Dune::FieldVector<ctype,dim>& testGrad, const Dune::FieldVector<ctype,dim>& ansatzGrad) const
{
Dune::FieldMatrix<ctype, dim, dim> res(0);
if (dim==2) {
res[0][1] = testGrad[0]*ansatzGrad[1]-ansatzGrad[0]*testGrad[1];
res[1][0] = -res[0][1];
......@@ -217,7 +217,7 @@ private:
}
//compute the cross product
Dune::FieldVector<ctype,dim> cross(0);
Dune::FieldVector<ctype,dim> cross(0);
for (int k=0; k<dim; k++)
cross[k]=testGrad[(k+1)%dim]*ansatzGrad[(k+2)%dim] -testGrad[(k+2)%dim]*ansatzGrad[(k+1)%dim];
......@@ -228,7 +228,7 @@ private:
res[i][j] = (cross*defGrad[k])*std::pow(-1,k);
res[j][i] = -res[i][j];
}
return res;
}
};
......
dune/elasticity/common/trustregionsolver.cc
View file @ 2fde2da4
......@@ -177,10 +177,10 @@ setup(const typename BasisType::GridView::Grid& grid,
typedef BasisType Basis;
bool isP1Basis = std::is_same<Basis,Dune::Functions::LagrangeBasis<typename Basis::GridView, 1> >::value;
using TransferOperatorType = typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType;
std::vector<std::shared_ptr<TruncatedCompressedMGTransfer<CorrectionType>>> transferOperators(isP1Basis ? numLevels-1 : numLevels);
// Here we set up the restriction of the leaf grid space into the leaf grid P1/Q1 space
if (not isP1Basis)
{
......@@ -318,7 +318,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
//Take the obstacles on the finest grid and give them to the multigrid solver, it will create a hierarchy for all coarser grids
trustRegionObstacles = trustRegion.obstacles();
#if ! HAVE_DUNE_PARMG
mgStep->setProblem(stiffnessMatrix, corr, rhs);
......
dune/elasticity/estimators/zienkiewiczzhuestimator.hh
View file @ 2fde2da4
......@@ -27,11 +27,11 @@ public:
: grid_(&grid), E_(E), nu_(nu)
{}
void estimate(const VectorType& x,
void estimate(const VectorType& x,
Dune::BlockVector<Dune::FieldVector<double,1> >& errorPerElement) {
using namespace Dune;
const auto& indexSet = grid_->leafIndexSet();
......@@ -46,8 +46,8 @@ public:
typedef BCRSMatrix<BlockType> MassMatrixType;
typedef P1NodalBasis<typename GridType::LeafGridView,double> P1Basis;
P1Basis p1Basis(grid_->leafGridView());
MassMatrixType massMatrix;
OperatorAssembler<P1Basis,P1Basis> globalAssembler(p1Basis,p1Basis);
......@@ -60,28 +60,28 @@ public:
#ifdef LUMP_MATRIX
/* BDMatrix<FieldMatrix<double,1,1> > invMassMatrix(baseSet.size());
invMassMatrix = 0;
for (int i=0; i<baseSet.size(); i++) {
for (int j=0; j<baseSet.size(); j++)
invMassMatrix[i][i] += massMatrix[i][j];
}
*/
BDMatrix<BlockType> invMassMatrix;
invMassMatrix = 0;
lumpMatrix(massMatrix,invMassMatrix);
invMassMatrix.invert();
#endif
#endif
BlockVector<SymmetricTensor<dim> > unscaledP1Stress(indexSet.size(dim));
unscaledP1Stress = 0;
typedef Dune::PQkLocalFiniteElementCache<typename GridType::ctype, double, dim, 1> FiniteElementCache;
FiniteElementCache cache;
FiniteElementCache cache;
for (const auto& e : elements(grid_->leafGridView())) {
// Element type
......@@ -110,10 +110,10 @@ public:
auto&& geometry = e.geometry();
for (size_t g=0; g<quad.size(); ++g) {
// pos of integration point
const FieldVector<double,dim>& quadPos = quad[g].position();
// eval jacobian inverse
auto&& jac = geometry.jacobianInverseTransposed(quadPos);
......@@ -124,38 +124,38 @@ public:
// Compute p0 stress at this quadrature node
// /////////////////////////////////////////////
SymmetricTensor<dim> p0Strain(0.0);
// evaluate gradients at Gauss points
std::vector<FieldMatrix<double,1,dim> > temp;
FieldVector<double,dim> grad;
finiteElement.localBasis().evaluateJacobian(quadPos,temp);
for (size_t i=0; i<finiteElement.localBasis().size(); i++) {
grad = 0;
jac.umv(temp[i][0],grad); // transform gradient to global coordinates
for (int j=0; j<dim; j++) {
// The deformation gradient of the shape function
FieldMatrix<double, dim, dim> Grad(0);
Grad[j] = grad;
/* Computes the linear strain tensor from the deformation gradient*/
SymmetricTensor<dim> scaledStrain;
computeStrain(Grad,scaledStrain);
scaledStrain *= localSolution[i][j];
p0Strain += scaledStrain;
}
}
// compute stress
SymmetricTensor<dim> p0Stress = hookeTimesStrain(p0Strain);
std::vector<FieldVector<double,1> >values;
finiteElement.localBasis().evaluateFunction(quadPos,values);
......@@ -163,15 +163,15 @@ public:
for (size_t row=0; row<finiteElement.localBasis().size(); row++) {
int globalRow = indexSet.subIndex(e,row,dim);
for (size_t rcomp=0; rcomp<p0Stress.size(); rcomp++) {
unscaledP1Stress[globalRow][rcomp] += p0Stress[rcomp]*values[row] * factor;
}
}
}
}
}
......@@ -182,7 +182,7 @@ public:
BlockVector<SymmetricTensor<dim> > elementP1Stress(indexSet.size(dim));
elementP1Stress = 0;
#ifdef LUMP_MATRIX
/*elementP1Stress = unscaledP1Stress;
for (int i=0; i<baseSet.size(); i++)
......@@ -192,20 +192,20 @@ public:
tmp = unscaledP1Stress;
invMassMatrix.mv(tmp,elementP1Stress);
#else
// Technicality: turn the matrix into a linear operator
MatrixAdapter<MassMatrixType, BlockVector<SymmetricTensor<dim> >, BlockVector<SymmetricTensor<dim> > > op(massMatrix);
// A preconditioner
SeqILU<MassMatrixType, BlockVector<SymmetricTensor<dim> >, BlockVector<SymmetricTensor<dim> > > ilu0(massMatrix,1.0);
// A preconditioned conjugate-gradient solver
int numIt = 100;
CGSolver<BlockVector<SymmetricTensor<dim> > > cg(op,ilu0,1E-16,numIt,0);
CGSolver<BlockVector<SymmetricTensor<dim> > > cg(op,ilu0,1E-16,numIt,0);
// Object storing some statistics about the solving process
InverseOperatorResult statistics;
// Solve!
cg.apply(elementP1Stress, unscaledP1Stress, statistics);
#endif
......@@ -216,7 +216,7 @@ public:
GeometryType gt = e.type();
// First order finite element
const typename FiniteElementCache::FiniteElementType& finiteElement = cache.get(gt);
const typename FiniteElementCache::FiniteElementType& finiteElement = cache.get(gt);
// /////////////////////////////////////////
// Extract the solution on this element
......@@ -238,22 +238,22 @@ public:
errorPerElement[indexSet.index(e)] = 0;
for (size_t g=0; g<quad.size(); ++g) {
// pos of integration point
const auto& quadPos = quad[g].position();
// eval jacobian inverse
const auto& jac = e.geometry().jacobianInverseTransposed(quadPos);
// weight of quadrature point
double weight = quad[g].weight();
// determinant of jacobian
double detjac = e.geometry().integrationElement(quadPos);
// Stresses at this quadrature point
SymmetricTensor<dim> p0Stress(0.0), p1Stress(0.0);
// evaluate gradients at Gauss points
std::vector<FieldMatrix<double,1,dim> >temp;
std::vector<FieldVector<double,1> >values;
......@@ -264,16 +264,16 @@ public:
// loop over test function number
for (size_t i=0; i<finiteElement.localBasis().size(); i++) {
grad = 0;
jac.umv(temp[i][0],grad); // transform gradient to global coordinates
for (int j=0; j<dim; j++) {
// Compute strain of shape functions
FieldMatrix<double, dim, dim> Grad(0);
Grad[j] = grad;
/* Computes the linear strain tensor from the deformation gradient*/
SymmetricTensor<dim> shapeFunctionStrain;
computeStrain(Grad,shapeFunctionStrain);
......@@ -282,7 +282,7 @@ public:
// Compute stress
p0Stress += hookeTimesStrain(shapeFunctionStrain);
}
}
SymmetricTensor<dim> scaledStress = elementP1Stress[i];
scaledStress *= values[i];
......@@ -297,7 +297,7 @@ public:
}
private:
void computeStrain(const Dune::FieldMatrix<double, dim, dim>& grad, SymmetricTensor<dim>& strain) const
{
for (int i=0; i<dim ; ++i)
......@@ -306,11 +306,11 @@ private:
for (int j=i+1; j<dim; ++j)
strain(i,j) = 0.5*(grad[i][j] + grad[j][i]);
}
}
SymmetricTensor<dim> hookeTimesStrain(const SymmetricTensor<dim>& strain) const {
SymmetricTensor<dim> stress = strain;
stress *= E_/(1+nu_);
......
dune/elasticity/materials/geomexactstvenantkirchhoffmaterial.hh
View file @ 2fde2da4
......@@ -14,7 +14,7 @@
/* \brief Class representing a hyperelastic geometrically exact St.Venant Kirchhoff material.
*
* \tparam Basis Global basis that is used for the spatial discretization.
* \tparam Basis Global basis that is used for the spatial discretization.
* (This is not nice but I need a LocalFiniteElement for the local Hessian assembler :-( )
*/
template <class Basis>
......@@ -111,14 +111,14 @@ public:
return 0.5*energy;
}
//! Return the local assembler of the first derivative of the strain energy
//! Return the local assembler of the first derivative of the strain energy
LinearisationAssembler& firstDerivative(std::shared_ptr<GridFunction> displace)
{
localLinearisation_.setConfiguration(displace);
return localLinearisation_;
}
//! Return the local assembler of the second derivative of the strain energy
//! Return the local assembler of the second derivative of the strain energy
HessianAssembler& secondDerivative(std::shared_ptr<GridFunction> displace)
{
localHessian_.setConfiguration(displace);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment