48 #include "CompressibleNonlinearElasticitySolver.hpp"
49 #include "LinearBasisFunction.hpp"
50 #include "QuadraticBasisFunction.hpp"
55 bool assembleJacobian)
58 assert(assembleResidual || assembleJacobian);
59 assert(this->mCurrentSolution.size()==this->mNumDofs);
73 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE> a_elem;
76 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE> a_elem_precond;
77 c_vector<double, STENCIL_SIZE> b_elem;
81 iter != this->mrQuadMesh.GetElementIteratorEnd();
88 #define COVERAGE_IGNORE
92 std::cout <<
"\r[" <<
PetscTools::GetMyRank() <<
"]: Element " << (*iter).GetIndex() <<
" of " << this->mrQuadMesh.GetNumElements() << std::flush;
94 #undef COVERAGE_IGNORE
96 AssembleOnElement(element, a_elem, a_elem_precond, b_elem, assembleResidual, assembleJacobian);
109 unsigned p_indices[STENCIL_SIZE];
110 for (
unsigned i=0; i<NUM_NODES_PER_ELEMENT; i++)
112 for (
unsigned j=0; j<DIM; j++)
118 if (assembleJacobian)
120 PetscMatTools::AddMultipleValues<STENCIL_SIZE>(this->mrJacobianMatrix, p_indices, a_elem);
121 PetscMatTools::AddMultipleValues<STENCIL_SIZE>(this->mPreconditionMatrix, p_indices, a_elem_precond);
124 if (assembleResidual)
126 PetscVecTools::AddMultipleValues<STENCIL_SIZE>(this->mResidualVector, p_indices, b_elem);
132 c_vector<double, BOUNDARY_STENCIL_SIZE> b_boundary_elem;
133 c_matrix<double, BOUNDARY_STENCIL_SIZE, BOUNDARY_STENCIL_SIZE> a_boundary_elem;
134 if (this->mrProblemDefinition.GetTractionBoundaryConditionType() != NO_TRACTIONS)
136 for (
unsigned bc_index=0; bc_index<this->mrProblemDefinition.rGetTractionBoundaryElements().size(); bc_index++)
138 BoundaryElement<DIM-1,DIM>& r_boundary_element = *(this->mrProblemDefinition.rGetTractionBoundaryElements()[bc_index]);
147 this->AssembleOnBoundaryElement(r_boundary_element, a_boundary_elem, b_boundary_elem, assembleResidual, assembleJacobian, bc_index);
149 unsigned p_indices[BOUNDARY_STENCIL_SIZE];
150 for (
unsigned i=0; i<NUM_NODES_PER_BOUNDARY_ELEMENT; i++)
152 for (
unsigned j=0; j<DIM; j++)
154 p_indices[DIM*i+j] = DIM*r_boundary_element.GetNodeGlobalIndex(i) + j;
158 if (assembleJacobian)
160 PetscMatTools::AddMultipleValues<BOUNDARY_STENCIL_SIZE>(this->mrJacobianMatrix, p_indices, a_boundary_elem);
161 PetscMatTools::AddMultipleValues<BOUNDARY_STENCIL_SIZE>(this->mPreconditionMatrix, p_indices, a_boundary_elem);
164 if (assembleResidual)
166 PetscVecTools::AddMultipleValues<BOUNDARY_STENCIL_SIZE>(this->mResidualVector, p_indices, b_boundary_elem);
171 this->FinishAssembleSystem(assembleResidual, assembleJacobian);
177 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE >& rAElem,
178 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE >& rAElemPrecond,
179 c_vector<double, STENCIL_SIZE>& rBElem,
180 bool assembleResidual,
181 bool assembleJacobian)
183 static c_matrix<double,DIM,DIM> jacobian;
184 static c_matrix<double,DIM,DIM> inverse_jacobian;
185 double jacobian_determinant;
187 this->mrQuadMesh.GetInverseJacobianForElement(rElement.
GetIndex(), jacobian, jacobian_determinant, inverse_jacobian);
189 if (assembleJacobian)
192 rAElemPrecond.clear();
195 if (assembleResidual)
201 static c_matrix<double,DIM,NUM_NODES_PER_ELEMENT> element_current_displacements;
202 for (
unsigned II=0; II<NUM_NODES_PER_ELEMENT; II++)
204 for (
unsigned JJ=0; JJ<DIM; JJ++)
206 element_current_displacements(JJ,II) = this->mCurrentSolution[DIM*rElement.
GetNodeGlobalIndex(II) + JJ];
211 static c_vector<double, NUM_VERTICES_PER_ELEMENT> linear_phi;
212 static c_vector<double, NUM_NODES_PER_ELEMENT> quad_phi;
213 static c_matrix<double, DIM, NUM_NODES_PER_ELEMENT> grad_quad_phi;
214 static c_matrix<double, NUM_NODES_PER_ELEMENT, DIM> trans_grad_quad_phi;
218 = this->mrProblemDefinition.GetCompressibleMaterialLaw(rElement.
GetIndex());
221 static c_matrix<double,DIM,DIM> grad_u;
223 static c_matrix<double,DIM,DIM> F;
224 static c_matrix<double,DIM,DIM> C;
225 static c_matrix<double,DIM,DIM> inv_C;
226 static c_matrix<double,DIM,DIM> inv_F;
227 static c_matrix<double,DIM,DIM> T;
229 static c_matrix<double,DIM,DIM> F_T;
230 static c_matrix<double,DIM,NUM_NODES_PER_ELEMENT> F_T_grad_quad_phi;
232 c_vector<double,DIM> body_force;
240 static c_matrix<double, DIM, NUM_NODES_PER_ELEMENT> temp_matrix;
241 static c_matrix<double,NUM_NODES_PER_ELEMENT,DIM> grad_quad_phi_times_invF;
243 if(this->mSetComputeAverageStressPerElement)
245 this->mAverageStressesPerElement[rElement.
GetIndex()] = zero_vector<double>(DIM*(DIM+1)/2);
249 for (
unsigned quadrature_index=0; quadrature_index < this->mpQuadratureRule->GetNumQuadPoints(); quadrature_index++)
252 unsigned current_quad_point_global_index = rElement.
GetIndex()*this->mpQuadratureRule->GetNumQuadPoints()
255 double wJ = jacobian_determinant * this->mpQuadratureRule->GetWeight(quadrature_index);
257 const ChastePoint<DIM>& quadrature_point = this->mpQuadratureRule->rGetQuadPoint(quadrature_index);
263 trans_grad_quad_phi = trans(grad_quad_phi);
266 if (assembleResidual)
268 switch (this->mrProblemDefinition.GetBodyForceType())
270 case FUNCTIONAL_BODY_FORCE:
272 c_vector<double,DIM> X = zero_vector<double>(DIM);
274 for (
unsigned node_index=0; node_index<NUM_VERTICES_PER_ELEMENT; node_index++)
276 X += linear_phi(node_index)*this->mrQuadMesh.GetNode( rElement.
GetNodeGlobalIndex(node_index) )->rGetLocation();
278 body_force = this->mrProblemDefinition.EvaluateBodyForceFunction(X, this->mCurrentTime);
281 case CONSTANT_BODY_FORCE:
283 body_force = this->mrProblemDefinition.GetConstantBodyForce();
292 grad_u = zero_matrix<double>(DIM,DIM);
293 for (
unsigned node_index=0; node_index<NUM_NODES_PER_ELEMENT; node_index++)
295 for (
unsigned i=0; i<DIM; i++)
297 for (
unsigned M=0; M<DIM; M++)
299 grad_u(i,M) += grad_quad_phi(M,node_index)*element_current_displacements(i,node_index);
305 for (
unsigned i=0; i<DIM; i++)
307 for (
unsigned M=0; M<DIM; M++)
309 F(i,M) = (i==M?1:0) + grad_u(i,M);
313 C = prod(trans(F),F);
318 this->SetupChangeOfBasisMatrix(rElement.
GetIndex(), current_quad_point_global_index);
322 if(this->mIncludeActiveTension)
326 this->AddActiveStressAndStressDerivative(C, rElement.
GetIndex(), current_quad_point_global_index,
327 T, dTdE, assembleJacobian);
330 if(this->mSetComputeAverageStressPerElement)
332 this->AddStressToAverageStressPerElement(T,rElement.
GetIndex());
336 if (assembleResidual)
339 F_T_grad_quad_phi = prod(F_T, grad_quad_phi);
341 for (
unsigned index=0; index<NUM_NODES_PER_ELEMENT*DIM; index++)
343 unsigned spatial_dim = index%DIM;
344 unsigned node_index = (index-spatial_dim)/DIM;
346 rBElem(index) += - this->mrProblemDefinition.GetDensity()
347 * body_force(spatial_dim)
348 * quad_phi(node_index)
352 rBElem(index) += F_T_grad_quad_phi(spatial_dim,node_index)
358 if (assembleJacobian)
361 grad_quad_phi_times_invF = prod(trans_grad_quad_phi, inv_F);
376 for (
unsigned M=0; M<DIM; M++)
378 for (
unsigned N=0; N<DIM; N++)
380 for (
unsigned P=0; P<DIM; P++)
382 for (
unsigned Q=0; Q<DIM; Q++)
385 dSdF(M,N,P,Q) = 0.5*(dTdE(M,N,P,Q) + dTdE(M,N,Q,P));
392 dTdE.template SetAsContractionOnSecondDimension<DIM>(F, dSdF);
395 dSdF.template SetAsContractionOnFourthDimension<DIM>(F, dTdE);
398 for (
unsigned M=0; M<DIM; M++)
400 for (
unsigned N=0; N<DIM; N++)
402 for (
unsigned i=0; i<DIM; i++)
404 dSdF(M,i,N,i) += T(M,N);
421 temp_tensor.template SetAsContractionOnFirstDimension<DIM>(trans_grad_quad_phi, dSdF);
422 dSdF_quad_quad.template SetAsContractionOnThirdDimension<DIM>(trans_grad_quad_phi, temp_tensor);
424 for (
unsigned index1=0; index1<NUM_NODES_PER_ELEMENT*DIM; index1++)
426 unsigned spatial_dim1 = index1%DIM;
427 unsigned node_index1 = (index1-spatial_dim1)/DIM;
429 for (
unsigned index2=0; index2<NUM_NODES_PER_ELEMENT*DIM; index2++)
431 unsigned spatial_dim2 = index2%DIM;
432 unsigned node_index2 = (index2-spatial_dim2)/DIM;
435 rAElem(index1,index2) += dSdF_quad_quad(node_index1,spatial_dim1,node_index2,spatial_dim2)
442 rAElemPrecond.clear();
443 if (assembleJacobian)
445 rAElemPrecond = rAElem;
448 if(this->mSetComputeAverageStressPerElement)
450 for(
unsigned i=0; i<DIM*(DIM+1)/2; i++)
452 this->mAverageStressesPerElement[rElement.
GetIndex()](i) /= this->mpQuadratureRule->GetNumQuadPoints();
462 std::string outputDirectory)
470 EXCEPTION(
"SolidMechanicsProblemDefinition object contains incompressible material laws");
ElementIterator GetElementIteratorBegin(bool skipDeletedElements=true)
boost::numeric::ublas::c_matrix< T, 1, 1 > Inverse(const boost::numeric::ublas::c_matrix< T, 1, 1 > &rM)
void AssembleSystem(bool assembleResidual, bool assembleJacobian)
unsigned GetNodeGlobalIndex(unsigned localIndex) const
#define EXCEPTION(message)
static void ComputeTransformedBasisFunctionDerivatives(const ChastePoint< ELEMENT_DIM > &rPoint, const c_matrix< double, ELEMENT_DIM, ELEMENT_DIM > &rInverseJacobian, c_matrix< double, ELEMENT_DIM,(ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2 > &rReturnValue)
bool OptionExists(const std::string &rOption)
bool GetOwnership() const
virtual ~CompressibleNonlinearElasticitySolver()
static void ComputeBasisFunctions(const ChastePoint< ELEMENT_DIM > &rPoint, c_vector< double,(ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2 > &rReturnValue)
CompressibilityType GetCompressibilityType()
static void ComputeBasisFunctions(const ChastePoint< ELEMENT_DIM > &rPoint, c_vector< double, ELEMENT_DIM+1 > &rReturnValue)
virtual void AssembleOnElement(Element< DIM, DIM > &rElement, c_matrix< double, STENCIL_SIZE, STENCIL_SIZE > &rAElem, c_matrix< double, STENCIL_SIZE, STENCIL_SIZE > &rAElemPrecond, c_vector< double, STENCIL_SIZE > &rBElem, bool assembleResidual, bool assembleJacobian)
unsigned GetIndex() const
static CommandLineArguments * Instance()
virtual void ComputeStressAndStressDerivative(c_matrix< double, DIM, DIM > &rC, c_matrix< double, DIM, DIM > &rInvC, double pressure, c_matrix< double, DIM, DIM > &rT, FourthOrderTensor< DIM, DIM, DIM, DIM > &rDTdE, bool computeDTdE)=0
void SetChangeOfBasisMatrix(c_matrix< double, DIM, DIM > &rChangeOfBasisMatrix)
CompressibleNonlinearElasticitySolver(AbstractTetrahedralMesh< DIM, DIM > &rQuadMesh, SolidMechanicsProblemDefinition< DIM > &rProblemDefinition, std::string outputDirectory)