diff --git a/dune/elasticity/assemblers/mooneyrivlinfunctionalassembler.hh b/dune/elasticity/assemblers/mooneyrivlinfunctionalassembler.hh index 8e0a87fb47c18270a6d660f251206a1e609d7c10..733b0da71811a98ba0dd02652a0ed5bec8e42d77 100644 --- a/dune/elasticity/assemblers/mooneyrivlinfunctionalassembler.hh +++ b/dune/elasticity/assemblers/mooneyrivlinfunctionalassembler.hh @@ -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 diff --git a/dune/elasticity/assemblers/mooneyrivlinoperatorassembler.hh b/dune/elasticity/assemblers/mooneyrivlinoperatorassembler.hh index 9e551fb4dcb462826cb3fced91c1a344215a9e5d..43f7eaf915125af32e4b685af1c9b3ec781fcba2 100644 --- a/dune/elasticity/assemblers/mooneyrivlinoperatorassembler.hh +++ b/dune/elasticity/assemblers/mooneyrivlinoperatorassembler.hh @@ -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; } }; diff --git a/dune/elasticity/assemblers/neohookefunctionalassembler.hh b/dune/elasticity/assemblers/neohookefunctionalassembler.hh index 6ec2f1a8834b1b95170ab5491f524ce4887c1873..fd12a21aef6d80183fb224aba4075f5ed39d0641 100644 --- a/dune/elasticity/assemblers/neohookefunctionalassembler.hh +++ b/dune/elasticity/assemblers/neohookefunctionalassembler.hh @@ -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]); } } diff --git a/dune/elasticity/assemblers/neohookeoperatorassembler.hh b/dune/elasticity/assemblers/neohookeoperatorassembler.hh index a34853d118c8b5a716bf4a7465d071eb00f63b37..0920e12f6400d3a2bea91a2f04f3366533f9edb4 100644 --- a/dune/elasticity/assemblers/neohookeoperatorassembler.hh +++ b/dune/elasticity/assemblers/neohookeoperatorassembler.hh @@ -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; } }; diff --git a/dune/elasticity/common/trustregionsolver.cc b/dune/elasticity/common/trustregionsolver.cc index 304e2841088887eb294ceae5fac6031b8b9f9bbc..432b73a029fe6abefed83cc548915a07a03d413f 100644 --- a/dune/elasticity/common/trustregionsolver.cc +++ b/dune/elasticity/common/trustregionsolver.cc @@ -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); diff --git a/dune/elasticity/estimators/zienkiewiczzhuestimator.hh b/dune/elasticity/estimators/zienkiewiczzhuestimator.hh index 38dd8e942dca508af32b7e480d37e0471282cff8..95d28438f360a06c12855def2c77d2b099bb9349 100644 --- a/dune/elasticity/estimators/zienkiewiczzhuestimator.hh +++ b/dune/elasticity/estimators/zienkiewiczzhuestimator.hh @@ -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_); diff --git a/dune/elasticity/materials/geomexactstvenantkirchhoffmaterial.hh b/dune/elasticity/materials/geomexactstvenantkirchhoffmaterial.hh index 84c0e879a4a1a6b13acabb214c6265fe32e514f6..f899884647e8b17ba007120adaa458c3ac302c31 100644 --- a/dune/elasticity/materials/geomexactstvenantkirchhoffmaterial.hh +++ b/dune/elasticity/materials/geomexactstvenantkirchhoffmaterial.hh @@ -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);