53 bool assembleJacobian)
56 assert(assembleResidual || assembleJacobian);
57 assert(this->mCurrentSolution.size()==this->mNumDofs);
71 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE> a_elem;
74 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE> a_elem_precond;
75 c_vector<double, STENCIL_SIZE> b_elem;
79 iter != this->mrQuadMesh.GetElementIteratorEnd();
86 std::cout <<
"\r[" <<
PetscTools::GetMyRank() <<
"]: Element " << (*iter).GetIndex() <<
" of " << this->mrQuadMesh.GetNumElements() << std::flush;
94 AssembleOnElement(element, a_elem, a_elem_precond, b_elem, assembleResidual, assembleJacobian);
113 unsigned p_indices[STENCIL_SIZE];
114 for (
unsigned i=0; i<NUM_NODES_PER_ELEMENT; i++)
116 for (
unsigned j=0; j<DIM; j++)
123 for (
unsigned i=0; i<NUM_VERTICES_PER_ELEMENT; i++)
129 p_indices[DIM*NUM_NODES_PER_ELEMENT + i] = (DIM+1)*vertex_index + DIM;
132 if (assembleJacobian)
134 PetscMatTools::AddMultipleValues<STENCIL_SIZE>(this->mrJacobianMatrix, p_indices, a_elem);
135 PetscMatTools::AddMultipleValues<STENCIL_SIZE>(this->mPreconditionMatrix, p_indices, a_elem_precond);
138 if (assembleResidual)
140 PetscVecTools::AddMultipleValues<STENCIL_SIZE>(this->mResidualVector, p_indices, b_elem);
146 c_vector<double, BOUNDARY_STENCIL_SIZE> b_boundary_elem;
147 c_matrix<double, BOUNDARY_STENCIL_SIZE, BOUNDARY_STENCIL_SIZE> a_boundary_elem;
149 if (this->mrProblemDefinition.GetTractionBoundaryConditionType() != NO_TRACTIONS)
151 for (
unsigned bc_index=0; bc_index<this->mrProblemDefinition.rGetTractionBoundaryElements().size(); bc_index++)
153 BoundaryElement<DIM-1,DIM>& r_boundary_element = *(this->mrProblemDefinition.rGetTractionBoundaryElements()[bc_index]);
161 this->AssembleOnBoundaryElement(r_boundary_element, a_boundary_elem, b_boundary_elem, assembleResidual, assembleJacobian, bc_index);
163 unsigned p_indices[BOUNDARY_STENCIL_SIZE];
164 for (
unsigned i=0; i<NUM_NODES_PER_BOUNDARY_ELEMENT; i++)
166 for (
unsigned j=0; j<DIM; j++)
169 p_indices[DIM*i+j] = (DIM+1)*r_boundary_element.GetNodeGlobalIndex(i) + j;
173 if (assembleJacobian)
175 PetscMatTools::AddMultipleValues<BOUNDARY_STENCIL_SIZE>(this->mrJacobianMatrix, p_indices, a_boundary_elem);
176 PetscMatTools::AddMultipleValues<BOUNDARY_STENCIL_SIZE>(this->mPreconditionMatrix, p_indices, a_boundary_elem);
179 if (assembleResidual)
181 PetscVecTools::AddMultipleValues<BOUNDARY_STENCIL_SIZE>(this->mResidualVector, p_indices, b_boundary_elem);
187 if (assembleResidual)
191 if (assembleJacobian)
197 if (assembleJacobian)
199 this->AddIdentityBlockForDummyPressureVariables(NONLINEAR_PROBLEM_APPLY_TO_EVERYTHING);
201 else if (assembleResidual)
203 this->AddIdentityBlockForDummyPressureVariables(NONLINEAR_PROBLEM_APPLY_TO_RESIDUAL_ONLY);
206 this->FinishAssembleSystem(assembleResidual, assembleJacobian);
212 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE >& rAElem,
213 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE >& rAElemPrecond,
214 c_vector<double, STENCIL_SIZE>& rBElem,
215 bool assembleResidual,
216 bool assembleJacobian)
218 static c_matrix<double,DIM,DIM> jacobian;
219 static c_matrix<double,DIM,DIM> inverse_jacobian;
220 double jacobian_determinant;
222 this->mrQuadMesh.GetInverseJacobianForElement(rElement.
GetIndex(), jacobian, jacobian_determinant, inverse_jacobian);
224 if (assembleJacobian)
227 rAElemPrecond.clear();
230 if (assembleResidual)
236 static c_matrix<double,DIM,NUM_NODES_PER_ELEMENT> element_current_displacements;
237 static c_vector<double,NUM_VERTICES_PER_ELEMENT> element_current_pressures;
238 for (
unsigned II=0; II<NUM_NODES_PER_ELEMENT; II++)
240 for (
unsigned JJ=0; JJ<DIM; JJ++)
243 element_current_displacements(JJ,II) = this->mCurrentSolution[(DIM+1)*rElement.
GetNodeGlobalIndex(II) + JJ];
248 for (
unsigned II=0; II<NUM_VERTICES_PER_ELEMENT; II++)
255 element_current_pressures(II) = this->mCurrentSolution[(DIM+1)*vertex_index + DIM];
259 static c_vector<double, NUM_VERTICES_PER_ELEMENT> linear_phi;
260 static c_vector<double, NUM_NODES_PER_ELEMENT> quad_phi;
261 static c_matrix<double, DIM, NUM_NODES_PER_ELEMENT> grad_quad_phi;
262 static c_matrix<double, NUM_NODES_PER_ELEMENT, DIM> trans_grad_quad_phi;
266 = this->mrProblemDefinition.GetIncompressibleMaterialLaw(rElement.
GetIndex());
268 static c_matrix<double,DIM,DIM> grad_u;
270 static c_matrix<double,DIM,DIM> F;
271 static c_matrix<double,DIM,DIM> C;
272 static c_matrix<double,DIM,DIM> inv_C;
273 static c_matrix<double,DIM,DIM> inv_F;
274 static c_matrix<double,DIM,DIM> T;
276 static c_matrix<double,DIM,DIM> F_T;
277 static c_matrix<double,DIM,NUM_NODES_PER_ELEMENT> F_T_grad_quad_phi;
279 c_vector<double,DIM> body_force;
287 static c_matrix<double, DIM, NUM_NODES_PER_ELEMENT> temp_matrix;
288 static c_matrix<double,NUM_NODES_PER_ELEMENT,DIM> grad_quad_phi_times_invF;
291 if (this->mSetComputeAverageStressPerElement)
293 this->mAverageStressesPerElement[rElement.
GetIndex()] = zero_vector<double>(DIM*(DIM+1)/2);
297 for (
unsigned quadrature_index=0; quadrature_index < this->mpQuadratureRule->GetNumQuadPoints(); quadrature_index++)
300 unsigned current_quad_point_global_index = rElement.
GetIndex()*this->mpQuadratureRule->GetNumQuadPoints()
303 double wJ = jacobian_determinant * this->mpQuadratureRule->GetWeight(quadrature_index);
305 const ChastePoint<DIM>& quadrature_point = this->mpQuadratureRule->rGetQuadPoint(quadrature_index);
311 trans_grad_quad_phi = trans(grad_quad_phi);
314 if (assembleResidual)
316 switch (this->mrProblemDefinition.GetBodyForceType())
318 case FUNCTIONAL_BODY_FORCE:
320 c_vector<double,DIM> X = zero_vector<double>(DIM);
322 for (
unsigned node_index=0; node_index<NUM_VERTICES_PER_ELEMENT; node_index++)
324 X += linear_phi(node_index)*this->mrQuadMesh.GetNode( rElement.
GetNodeGlobalIndex(node_index) )->rGetLocation();
326 body_force = this->mrProblemDefinition.EvaluateBodyForceFunction(X, this->mCurrentTime);
329 case CONSTANT_BODY_FORCE:
331 body_force = this->mrProblemDefinition.GetConstantBodyForce();
340 grad_u = zero_matrix<double>(DIM,DIM);
342 for (
unsigned node_index=0; node_index<NUM_NODES_PER_ELEMENT; node_index++)
344 for (
unsigned i=0; i<DIM; i++)
346 for (
unsigned M=0; M<DIM; M++)
348 grad_u(i,M) += grad_quad_phi(M,node_index)*element_current_displacements(i,node_index);
354 for (
unsigned vertex_index=0; vertex_index<NUM_VERTICES_PER_ELEMENT; vertex_index++)
356 pressure += linear_phi(vertex_index)*element_current_pressures(vertex_index);
360 for (
unsigned i=0; i<DIM; i++)
362 for (
unsigned M=0; M<DIM; M++)
364 F(i,M) = (i==M?1:0) + grad_u(i,M);
368 C = prod(trans(F),F);
375 this->SetupChangeOfBasisMatrix(rElement.
GetIndex(), current_quad_point_global_index);
379 if (this->mIncludeActiveTension)
383 this->AddActiveStressAndStressDerivative(C, rElement.
GetIndex(), current_quad_point_global_index,
384 T, dTdE, assembleJacobian);
387 if (this->mSetComputeAverageStressPerElement)
389 this->AddStressToAverageStressPerElement(T,rElement.
GetIndex());
393 if (assembleResidual)
396 F_T_grad_quad_phi = prod(F_T, grad_quad_phi);
398 for (
unsigned index=0; index<NUM_NODES_PER_ELEMENT*DIM; index++)
400 unsigned spatial_dim = index%DIM;
401 unsigned node_index = (index-spatial_dim)/DIM;
403 rBElem(index) += - this->mrProblemDefinition.GetDensity()
404 * body_force(spatial_dim)
405 * quad_phi(node_index)
409 rBElem(index) += F_T_grad_quad_phi(spatial_dim,node_index)
413 for (
unsigned vertex_index=0; vertex_index<NUM_VERTICES_PER_ELEMENT; vertex_index++)
415 rBElem( NUM_NODES_PER_ELEMENT*DIM + vertex_index ) += linear_phi(vertex_index)
422 if (assembleJacobian)
425 grad_quad_phi_times_invF = prod(trans_grad_quad_phi, inv_F);
440 for (
unsigned M=0; M<DIM; M++)
442 for (
unsigned N=0; N<DIM; N++)
444 for (
unsigned P=0; P<DIM; P++)
446 for (
unsigned Q=0; Q<DIM; Q++)
449 dSdF(M,N,P,Q) = 0.5*(dTdE(M,N,P,Q) + dTdE(M,N,Q,P));
456 dTdE.template SetAsContractionOnSecondDimension<DIM>(F, dSdF);
459 dSdF.template SetAsContractionOnFourthDimension<DIM>(F, dTdE);
462 for (
unsigned M=0; M<DIM; M++)
464 for (
unsigned N=0; N<DIM; N++)
466 for (
unsigned i=0; i<DIM; i++)
468 dSdF(M,i,N,i) += T(M,N);
485 temp_tensor.template SetAsContractionOnFirstDimension<DIM>( trans_grad_quad_phi, dSdF );
486 dSdF_quad_quad.template SetAsContractionOnThirdDimension<DIM>( trans_grad_quad_phi, temp_tensor );
488 for (
unsigned index1=0; index1<NUM_NODES_PER_ELEMENT*DIM; index1++)
490 unsigned spatial_dim1 = index1%DIM;
491 unsigned node_index1 = (index1-spatial_dim1)/DIM;
494 for (
unsigned index2=0; index2<NUM_NODES_PER_ELEMENT*DIM; index2++)
496 unsigned spatial_dim2 = index2%DIM;
497 unsigned node_index2 = (index2-spatial_dim2)/DIM;
500 rAElem(index1,index2) += dSdF_quad_quad(node_index1,spatial_dim1,node_index2,spatial_dim2)
504 for (
unsigned vertex_index=0; vertex_index<NUM_VERTICES_PER_ELEMENT; vertex_index++)
506 unsigned index2 = NUM_NODES_PER_ELEMENT*DIM + vertex_index;
509 rAElem(index1,index2) += - grad_quad_phi_times_invF(node_index1,spatial_dim1)
510 * linear_phi(vertex_index)
515 for (
unsigned vertex_index=0; vertex_index<NUM_VERTICES_PER_ELEMENT; vertex_index++)
517 unsigned index1 = NUM_NODES_PER_ELEMENT*DIM + vertex_index;
519 for (
unsigned index2=0; index2<NUM_NODES_PER_ELEMENT*DIM; index2++)
521 unsigned spatial_dim2 = index2%DIM;
522 unsigned node_index2 = (index2-spatial_dim2)/DIM;
525 rAElem(index1,index2) += detF
526 * grad_quad_phi_times_invF(node_index2,spatial_dim2)
527 * linear_phi(vertex_index)
537 for (
unsigned vertex_index2=0; vertex_index2<NUM_VERTICES_PER_ELEMENT; vertex_index2++)
539 unsigned index2 = NUM_NODES_PER_ELEMENT*DIM + vertex_index2;
540 rAElemPrecond(index1,index2) += linear_phi(vertex_index)
541 * linear_phi(vertex_index2)
548 if (assembleJacobian)
550 if (this->mPetscDirectSolve)
556 rAElemPrecond = rAElemPrecond + rAElem;
567 rAElemPrecond = rAElemPrecond + rAElem;
569 for (
unsigned i=NUM_NODES_PER_ELEMENT*DIM; i<STENCIL_SIZE; i++)
571 for (
unsigned j=0; j<NUM_NODES_PER_ELEMENT*DIM; j++)
573 rAElemPrecond(i,j) = 0.0;
579 if (this->mSetComputeAverageStressPerElement)
581 for (
unsigned i=0; i<DIM*(DIM+1)/2; i++)
583 this->mAverageStressesPerElement[rElement.
GetIndex()](i) /= this->mpQuadratureRule->GetNumQuadPoints();