00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 #include "IncompressibleNonlinearElasticitySolver.hpp"
00040 #include "LinearBasisFunction.hpp"
00041 #include "QuadraticBasisFunction.hpp"
00042 #include <algorithm>
00043
00044 template<size_t DIM>
00045 void IncompressibleNonlinearElasticitySolver<DIM>::AssembleSystem(bool assembleResidual,
00046 bool assembleJacobian)
00047 {
00048
00049 assert(assembleResidual || assembleJacobian);
00050 assert(this->mCurrentSolution.size()==this->mNumDofs);
00051
00052
00053 if (assembleResidual)
00054 {
00055 PetscVecTools::Zero(this->mResidualVector);
00056 }
00057 if (assembleJacobian)
00058 {
00059 PetscMatTools::Zero(this->mJacobianMatrix);
00060 PetscMatTools::Zero(this->mPreconditionMatrix);
00061 }
00062
00063 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE> a_elem;
00064
00065
00066 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE> a_elem_precond;
00067 c_vector<double, STENCIL_SIZE> b_elem;
00068
00069
00070 for (typename AbstractTetrahedralMesh<DIM, DIM>::ElementIterator iter = this->mrQuadMesh.GetElementIteratorBegin();
00071 iter != this->mrQuadMesh.GetElementIteratorEnd();
00072 ++iter)
00073 {
00074 #ifdef MECH_VERY_VERBOSE
00075 if (assembleJacobian)
00076 {
00077 std::cout << "\r[" << PetscTools::GetMyRank() << "]: Element " << (*iter).GetIndex() << " of " << this->mrQuadMesh.GetNumElements() << std::flush;
00078 }
00079 #endif
00080
00081 Element<DIM, DIM>& element = *iter;
00082
00083 if (element.GetOwnership() == true)
00084 {
00085 AssembleOnElement(element, a_elem, a_elem_precond, b_elem, assembleResidual, assembleJacobian);
00086
00090
00091
00092
00093
00094
00095
00096
00097
00098 unsigned p_indices[STENCIL_SIZE];
00099 for (unsigned i=0; i<NUM_NODES_PER_ELEMENT; i++)
00100 {
00101 for (unsigned j=0; j<DIM; j++)
00102 {
00103 p_indices[DIM*i+j] = DIM*element.GetNodeGlobalIndex(i) + j;
00104 }
00105 }
00106
00107 for (unsigned i=0; i<NUM_VERTICES_PER_ELEMENT; i++)
00108 {
00109
00110
00111 unsigned vertex_index = element.GetNodeGlobalIndex(i);
00112 assert(vertex_index < this->mrQuadMesh.GetNumVertices());
00113
00114
00115
00116
00117
00118 p_indices[DIM*NUM_NODES_PER_ELEMENT + i] = DIM*this->mrQuadMesh.GetNumNodes() + vertex_index;
00119 }
00120
00121 if (assembleJacobian)
00122 {
00123 PetscMatTools::AddMultipleValues<STENCIL_SIZE>(this->mJacobianMatrix, p_indices, a_elem);
00124 PetscMatTools::AddMultipleValues<STENCIL_SIZE>(this->mPreconditionMatrix, p_indices, a_elem_precond);
00125 }
00126
00127 if (assembleResidual)
00128 {
00129 PetscVecTools::AddMultipleValues<STENCIL_SIZE>(this->mResidualVector, p_indices, b_elem);
00130 }
00131 }
00132 }
00133
00134
00135 c_vector<double, BOUNDARY_STENCIL_SIZE> b_boundary_elem;
00136 c_matrix<double, BOUNDARY_STENCIL_SIZE, BOUNDARY_STENCIL_SIZE> a_boundary_elem;
00137 if (this->mrProblemDefinition.GetTractionBoundaryConditionType() != NO_TRACTIONS)
00138 {
00139 for (unsigned bc_index=0; bc_index<this->mrProblemDefinition.rGetTractionBoundaryElements().size(); bc_index++)
00140 {
00141 BoundaryElement<DIM-1,DIM>& r_boundary_element = *(this->mrProblemDefinition.rGetTractionBoundaryElements()[bc_index]);
00142 AssembleOnBoundaryElement(r_boundary_element, a_boundary_elem, b_boundary_elem, assembleResidual, assembleJacobian, bc_index);
00143
00144 unsigned p_indices[BOUNDARY_STENCIL_SIZE];
00145 for (unsigned i=0; i<NUM_NODES_PER_BOUNDARY_ELEMENT; i++)
00146 {
00147 for (unsigned j=0; j<DIM; j++)
00148 {
00149 p_indices[DIM*i+j] = DIM*r_boundary_element.GetNodeGlobalIndex(i) + j;
00150 }
00151 }
00152
00153 for (unsigned i=0; i<DIM ; i++)
00154 {
00155 p_indices[DIM*NUM_NODES_PER_BOUNDARY_ELEMENT + i] = DIM*this->mrQuadMesh.GetNumNodes() + r_boundary_element.GetNodeGlobalIndex(i);
00156 }
00157
00158 if (assembleJacobian)
00159 {
00160 PetscMatTools::AddMultipleValues<BOUNDARY_STENCIL_SIZE>(this->mJacobianMatrix, p_indices, a_boundary_elem);
00161 PetscMatTools::AddMultipleValues<BOUNDARY_STENCIL_SIZE>(this->mPreconditionMatrix, p_indices, a_boundary_elem);
00162 }
00163
00164 if (assembleResidual)
00165 {
00166 PetscVecTools::AddMultipleValues<BOUNDARY_STENCIL_SIZE>(this->mResidualVector, p_indices, b_boundary_elem);
00167 }
00168
00169
00170 if (DIM==2)
00171 {
00172 assert(8==BOUNDARY_STENCIL_SIZE);
00173 assert(b_boundary_elem(6)==0);
00174 assert(b_boundary_elem(7)==0);
00175 }
00176 }
00177 }
00178
00179 this->FinishAssembleSystem(assembleResidual, assembleJacobian);
00180 }
00181
00182 template<size_t DIM>
00183 void IncompressibleNonlinearElasticitySolver<DIM>::AssembleOnElement(
00184 Element<DIM, DIM>& rElement,
00185 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE >& rAElem,
00186 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE >& rAElemPrecond,
00187 c_vector<double, STENCIL_SIZE>& rBElem,
00188 bool assembleResidual,
00189 bool assembleJacobian)
00190 {
00191 static c_matrix<double,DIM,DIM> jacobian;
00192 static c_matrix<double,DIM,DIM> inverse_jacobian;
00193 double jacobian_determinant;
00194
00195 this->mrQuadMesh.GetInverseJacobianForElement(rElement.GetIndex(), jacobian, jacobian_determinant, inverse_jacobian);
00196
00197 if (assembleJacobian)
00198 {
00199 rAElem.clear();
00200 rAElemPrecond.clear();
00201 }
00202
00203 if (assembleResidual)
00204 {
00205 rBElem.clear();
00206 }
00207
00208
00209 static c_matrix<double,DIM,NUM_NODES_PER_ELEMENT> element_current_displacements;
00210 static c_vector<double,NUM_VERTICES_PER_ELEMENT> element_current_pressures;
00211 for (unsigned II=0; II<NUM_NODES_PER_ELEMENT; II++)
00212 {
00213 for (unsigned JJ=0; JJ<DIM; JJ++)
00214 {
00215 element_current_displacements(JJ,II) = this->mCurrentSolution[DIM*rElement.GetNodeGlobalIndex(II) + JJ];
00216 }
00217 }
00218
00219
00220 for (unsigned II=0; II<NUM_VERTICES_PER_ELEMENT; II++)
00221 {
00222
00223
00224 unsigned vertex_index = rElement.GetNodeGlobalIndex(II);
00225 assert(vertex_index < this->mrQuadMesh.GetNumVertices());
00226
00227
00228
00229
00230
00231 element_current_pressures(II) = this->mCurrentSolution[DIM*this->mrQuadMesh.GetNumNodes() + vertex_index];
00232 }
00233
00234
00235 static c_vector<double, NUM_VERTICES_PER_ELEMENT> linear_phi;
00236 static c_vector<double, NUM_NODES_PER_ELEMENT> quad_phi;
00237 static c_matrix<double, DIM, NUM_NODES_PER_ELEMENT> grad_quad_phi;
00238 static c_matrix<double, NUM_NODES_PER_ELEMENT, DIM> trans_grad_quad_phi;
00239
00240
00241 AbstractIncompressibleMaterialLaw<DIM>* p_material_law
00242 = this->mrProblemDefinition.GetIncompressibleMaterialLaw(rElement.GetIndex());
00243
00244 static c_matrix<double,DIM,DIM> grad_u;
00245
00246 static c_matrix<double,DIM,DIM> F;
00247 static c_matrix<double,DIM,DIM> C;
00248 static c_matrix<double,DIM,DIM> inv_C;
00249 static c_matrix<double,DIM,DIM> inv_F;
00250 static c_matrix<double,DIM,DIM> T;
00251
00252 static c_matrix<double,DIM,DIM> F_T;
00253 static c_matrix<double,DIM,NUM_NODES_PER_ELEMENT> F_T_grad_quad_phi;
00254
00255 c_vector<double,DIM> body_force;
00256
00257 static FourthOrderTensor<DIM,DIM,DIM,DIM> dTdE;
00258 static FourthOrderTensor<DIM,DIM,DIM,DIM> dSdF;
00259
00260 static FourthOrderTensor<NUM_NODES_PER_ELEMENT,DIM,DIM,DIM> temp_tensor;
00261 static FourthOrderTensor<NUM_NODES_PER_ELEMENT,DIM,NUM_NODES_PER_ELEMENT,DIM> dSdF_quad_quad;
00262
00263 static c_matrix<double, DIM, NUM_NODES_PER_ELEMENT> temp_matrix;
00264 static c_matrix<double,NUM_NODES_PER_ELEMENT,DIM> grad_quad_phi_times_invF;
00265
00266
00267 for (unsigned quadrature_index=0; quadrature_index < this->mpQuadratureRule->GetNumQuadPoints(); quadrature_index++)
00268 {
00269
00270 unsigned current_quad_point_global_index = rElement.GetIndex()*this->mpQuadratureRule->GetNumQuadPoints()
00271 + quadrature_index;
00272
00273 double wJ = jacobian_determinant * this->mpQuadratureRule->GetWeight(quadrature_index);
00274
00275 const ChastePoint<DIM>& quadrature_point = this->mpQuadratureRule->rGetQuadPoint(quadrature_index);
00276
00277
00278 LinearBasisFunction<DIM>::ComputeBasisFunctions(quadrature_point, linear_phi);
00279 QuadraticBasisFunction<DIM>::ComputeBasisFunctions(quadrature_point, quad_phi);
00280 QuadraticBasisFunction<DIM>::ComputeTransformedBasisFunctionDerivatives(quadrature_point, inverse_jacobian, grad_quad_phi);
00281 trans_grad_quad_phi = trans(grad_quad_phi);
00282
00283
00284 if (assembleResidual)
00285 {
00286 switch (this->mrProblemDefinition.GetBodyForceType())
00287 {
00288 case FUNCTIONAL_BODY_FORCE:
00289 {
00290 c_vector<double,DIM> X = zero_vector<double>(DIM);
00291
00292 for (unsigned node_index=0; node_index<NUM_VERTICES_PER_ELEMENT; node_index++)
00293 {
00294 X += linear_phi(node_index)*this->mrQuadMesh.GetNode( rElement.GetNodeGlobalIndex(node_index) )->rGetLocation();
00295 }
00296 body_force = this->mrProblemDefinition.EvaluateBodyForceFunction(X, this->mCurrentTime);
00297 break;
00298 }
00299 case CONSTANT_BODY_FORCE:
00300 {
00301 body_force = this->mrProblemDefinition.GetConstantBodyForce();
00302 break;
00303 }
00304 default:
00305 NEVER_REACHED;
00306 }
00307 }
00308
00309
00310 grad_u = zero_matrix<double>(DIM,DIM);
00311
00312 for (unsigned node_index=0; node_index<NUM_NODES_PER_ELEMENT; node_index++)
00313 {
00314 for (unsigned i=0; i<DIM; i++)
00315 {
00316 for (unsigned M=0; M<DIM; M++)
00317 {
00318 grad_u(i,M) += grad_quad_phi(M,node_index)*element_current_displacements(i,node_index);
00319 }
00320 }
00321 }
00322
00323 double pressure = 0;
00324 for (unsigned vertex_index=0; vertex_index<NUM_VERTICES_PER_ELEMENT; vertex_index++)
00325 {
00326 pressure += linear_phi(vertex_index)*element_current_pressures(vertex_index);
00327 }
00328
00329
00330 for (unsigned i=0; i<DIM; i++)
00331 {
00332 for (unsigned M=0; M<DIM; M++)
00333 {
00334 F(i,M) = (i==M?1:0) + grad_u(i,M);
00335 }
00336 }
00337
00338 C = prod(trans(F),F);
00339 inv_C = Inverse(C);
00340 inv_F = Inverse(F);
00341
00342 double detF = Determinant(F);
00343
00344 this->ComputeStressAndStressDerivative(p_material_law, C, inv_C, pressure, rElement.GetIndex(), current_quad_point_global_index,
00345 T, dTdE, assembleJacobian);
00346
00347
00348 if (assembleResidual)
00349 {
00350 F_T = prod(F,T);
00351 F_T_grad_quad_phi = prod(F_T, grad_quad_phi);
00352
00353 for (unsigned index=0; index<NUM_NODES_PER_ELEMENT*DIM; index++)
00354 {
00355 unsigned spatial_dim = index%DIM;
00356 unsigned node_index = (index-spatial_dim)/DIM;
00357
00358 rBElem(index) += - this->mrProblemDefinition.GetDensity()
00359 * body_force(spatial_dim)
00360 * quad_phi(node_index)
00361 * wJ;
00362
00363
00364 rBElem(index) += F_T_grad_quad_phi(spatial_dim,node_index)
00365 * wJ;
00366 }
00367
00368 for (unsigned vertex_index=0; vertex_index<NUM_VERTICES_PER_ELEMENT; vertex_index++)
00369 {
00370 rBElem( NUM_NODES_PER_ELEMENT*DIM + vertex_index ) += linear_phi(vertex_index)
00371 * (detF - 1)
00372 * wJ;
00373 }
00374 }
00375
00376
00377 if (assembleJacobian)
00378 {
00379
00380 grad_quad_phi_times_invF = prod(trans_grad_quad_phi, inv_F);
00381
00383
00384
00385
00386
00387
00388
00389
00390
00391
00393
00394
00395 for (unsigned M=0; M<DIM; M++)
00396 {
00397 for (unsigned N=0; N<DIM; N++)
00398 {
00399 for (unsigned P=0; P<DIM; P++)
00400 {
00401 for (unsigned Q=0; Q<DIM; Q++)
00402 {
00403
00404 dSdF(M,N,P,Q) = 0.5*(dTdE(M,N,P,Q) + dTdE(M,N,Q,P));
00405 }
00406 }
00407 }
00408 }
00409
00410
00411 dTdE.template SetAsContractionOnSecondDimension<DIM>(F, dSdF);
00412
00413
00414 dSdF.template SetAsContractionOnFourthDimension<DIM>(F, dTdE);
00415
00416
00417 for (unsigned M=0; M<DIM; M++)
00418 {
00419 for (unsigned N=0; N<DIM; N++)
00420 {
00421 for (unsigned i=0; i<DIM; i++)
00422 {
00423 dSdF(M,i,N,i) += T(M,N);
00424 }
00425 }
00426 }
00427
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00440 temp_tensor.template SetAsContractionOnFirstDimension<DIM>( trans_grad_quad_phi, dSdF );
00441 dSdF_quad_quad.template SetAsContractionOnThirdDimension<DIM>( trans_grad_quad_phi, temp_tensor );
00442
00443 for (unsigned index1=0; index1<NUM_NODES_PER_ELEMENT*DIM; index1++)
00444 {
00445 unsigned spatial_dim1 = index1%DIM;
00446 unsigned node_index1 = (index1-spatial_dim1)/DIM;
00447
00448
00449 for (unsigned index2=0; index2<NUM_NODES_PER_ELEMENT*DIM; index2++)
00450 {
00451 unsigned spatial_dim2 = index2%DIM;
00452 unsigned node_index2 = (index2-spatial_dim2)/DIM;
00453
00454
00455 rAElem(index1,index2) += dSdF_quad_quad(node_index1,spatial_dim1,node_index2,spatial_dim2)
00456 * wJ;
00457 }
00458
00459 for (unsigned vertex_index=0; vertex_index<NUM_VERTICES_PER_ELEMENT; vertex_index++)
00460 {
00461 unsigned index2 = NUM_NODES_PER_ELEMENT*DIM + vertex_index;
00462
00463
00464 rAElem(index1,index2) += - grad_quad_phi_times_invF(node_index1,spatial_dim1)
00465 * linear_phi(vertex_index)
00466 * wJ;
00467 }
00468 }
00469
00470 for (unsigned vertex_index=0; vertex_index<NUM_VERTICES_PER_ELEMENT; vertex_index++)
00471 {
00472 unsigned index1 = NUM_NODES_PER_ELEMENT*DIM + vertex_index;
00473
00474 for (unsigned index2=0; index2<NUM_NODES_PER_ELEMENT*DIM; index2++)
00475 {
00476 unsigned spatial_dim2 = index2%DIM;
00477 unsigned node_index2 = (index2-spatial_dim2)/DIM;
00478
00479
00480 rAElem(index1,index2) += detF
00481 * grad_quad_phi_times_invF(node_index2,spatial_dim2)
00482 * linear_phi(vertex_index)
00483 * wJ;
00484 }
00485
00487
00488
00489
00490
00492 for (unsigned vertex_index2=0; vertex_index2<NUM_VERTICES_PER_ELEMENT; vertex_index2++)
00493 {
00494 unsigned index2 = NUM_NODES_PER_ELEMENT*DIM + vertex_index2;
00495 rAElemPrecond(index1,index2) += linear_phi(vertex_index)
00496 * linear_phi(vertex_index2)
00497 * wJ;
00498 }
00499 }
00500 }
00501 }
00502
00503 if (assembleJacobian)
00504 {
00505
00506
00507
00508
00509
00510
00511
00512 rAElemPrecond = rAElemPrecond + rAElem;
00513 for (unsigned i=NUM_NODES_PER_ELEMENT*DIM; i<STENCIL_SIZE; i++)
00514 {
00515 for (unsigned j=0; j<NUM_NODES_PER_ELEMENT*DIM; j++)
00516 {
00517 rAElemPrecond(i,j) = 0.0;
00518 }
00519 }
00520 }
00521 }
00522
00523 template<size_t DIM>
00524 void IncompressibleNonlinearElasticitySolver<DIM>::AssembleOnBoundaryElement(
00525 BoundaryElement<DIM-1,DIM>& rBoundaryElement,
00526 c_matrix<double,BOUNDARY_STENCIL_SIZE,BOUNDARY_STENCIL_SIZE>& rAelem,
00527 c_vector<double,BOUNDARY_STENCIL_SIZE>& rBelem,
00528 bool assembleResidual,
00529 bool assembleJacobian,
00530 unsigned boundaryConditionIndex)
00531 {
00532 rAelem.clear();
00533 rBelem.clear();
00534
00535 if (assembleJacobian && !assembleResidual)
00536 {
00537
00538 return;
00539 }
00540
00541 c_vector<double, DIM> weighted_direction;
00542 double jacobian_determinant;
00543
00544 this->mrQuadMesh.GetWeightedDirectionForBoundaryElement(rBoundaryElement.GetIndex(), weighted_direction, jacobian_determinant);
00545
00546
00547
00548
00549
00550
00551 c_vector<double,DIM> deformed_normal;
00552 if (this->mrProblemDefinition.GetTractionBoundaryConditionType()==PRESSURE_ON_DEFORMED)
00553 {
00554
00555 static std::vector<c_vector<double,DIM> > element_current_displacements(DIM);
00556 for (unsigned II=0; II<DIM; II++)
00557 {
00558 for (unsigned JJ=0; JJ<DIM; JJ++)
00559 {
00560 element_current_displacements[II](JJ) = this->mCurrentSolution[DIM*rBoundaryElement.GetNodeGlobalIndex(II) + JJ];
00561 }
00562 }
00563
00564
00565 this->mDeformedBoundaryElement.ApplyUndeformedElementAndDisplacement(&rBoundaryElement, element_current_displacements);
00566
00567 this->mDeformedBoundaryElement.CalculateWeightedDirection(weighted_direction, jacobian_determinant);
00568
00569 deformed_normal = this->mDeformedBoundaryElement.CalculateNormal();
00570 }
00571
00572
00573
00574 c_vector<double,NUM_NODES_PER_BOUNDARY_ELEMENT> phi;
00575
00576 for (unsigned quad_index=0; quad_index<this->mpBoundaryQuadratureRule->GetNumQuadPoints(); quad_index++)
00577 {
00578 double wJ = jacobian_determinant * this->mpBoundaryQuadratureRule->GetWeight(quad_index);
00579
00580 const ChastePoint<DIM-1>& quad_point = this->mpBoundaryQuadratureRule->rGetQuadPoint(quad_index);
00581
00582 QuadraticBasisFunction<DIM-1>::ComputeBasisFunctions(quad_point, phi);
00583
00584
00585
00586 c_vector<double,DIM> traction = zero_vector<double>(DIM);
00587
00588 switch (this->mrProblemDefinition.GetTractionBoundaryConditionType())
00589 {
00590 case FUNCTIONAL_TRACTION:
00591 {
00592 c_vector<double,DIM> X = zero_vector<double>(DIM);
00593 for (unsigned node_index=0; node_index<NUM_NODES_PER_BOUNDARY_ELEMENT; node_index++)
00594 {
00595 X += phi(node_index)*this->mrQuadMesh.GetNode( rBoundaryElement.GetNodeGlobalIndex(node_index) )->rGetLocation();
00596 }
00597 traction = this->mrProblemDefinition.EvaluateTractionFunction(X, this->mCurrentTime);
00598 break;
00599 }
00600 case ELEMENTWISE_TRACTION:
00601 {
00602 traction = this->mrProblemDefinition.rGetElementwiseTractions()[boundaryConditionIndex];
00603 break;
00604 }
00605 case PRESSURE_ON_DEFORMED:
00606 {
00607
00608 traction = this->mrProblemDefinition.GetNormalPressure()*deformed_normal;
00609 break;
00610 }
00611 default:
00612 NEVER_REACHED;
00613 }
00614
00615
00616 for (unsigned index=0; index<NUM_NODES_PER_BOUNDARY_ELEMENT*DIM; index++)
00617 {
00618 unsigned spatial_dim = index%DIM;
00619 unsigned node_index = (index-spatial_dim)/DIM;
00620
00621 assert(node_index < NUM_NODES_PER_BOUNDARY_ELEMENT);
00622
00623 rBelem(index) -= traction(spatial_dim)
00624 * phi(node_index)
00625 * wJ;
00626 }
00627 }
00628 }
00629
00630 template<size_t DIM>
00631 void IncompressibleNonlinearElasticitySolver<DIM>::FormInitialGuess()
00632 {
00633 this->mCurrentSolution.resize(this->mNumDofs, 0.0);
00634
00635 for (unsigned i=0; i<this->mrQuadMesh.GetNumElements(); i++)
00636 {
00637 double zero_strain_pressure
00638 = this->mrProblemDefinition.GetIncompressibleMaterialLaw(i)->GetZeroStrainPressure();
00639
00640
00641
00642 for (unsigned j=0; j<NUM_VERTICES_PER_ELEMENT; j++)
00643 {
00644
00645
00646 unsigned vertex_index = this->mrQuadMesh.GetElement(i)->GetNodeGlobalIndex(j);
00647 assert(vertex_index < this->mrQuadMesh.GetNumVertices());
00648
00649
00650
00651
00652
00653 this->mCurrentSolution[ DIM*this->mrQuadMesh.GetNumNodes() + vertex_index ] = zero_strain_pressure;
00654 }
00655 }
00656 }
00657
00658 template<size_t DIM>
00659 IncompressibleNonlinearElasticitySolver<DIM>::IncompressibleNonlinearElasticitySolver(
00660 QuadraticMesh<DIM>& rQuadMesh,
00661 SolidMechanicsProblemDefinition<DIM>& rProblemDefinition,
00662 std::string outputDirectory)
00663 : AbstractNonlinearElasticitySolver<DIM>(rQuadMesh,
00664 rProblemDefinition,
00665 outputDirectory,
00666 INCOMPRESSIBLE)
00667 {
00668 if(rProblemDefinition.GetCompressibilityType() != INCOMPRESSIBLE)
00669 {
00670 EXCEPTION("SolidMechanicsProblemDefinition object contains compressible material laws");
00671 }
00672
00673 this->Initialise();
00674 FormInitialGuess();
00675 }
00676
00677
00678 template<size_t DIM>
00679 std::vector<double>& IncompressibleNonlinearElasticitySolver<DIM>::rGetPressures()
00680 {
00681 mPressures.clear();
00682 mPressures.resize(this->mrQuadMesh.GetNumVertices());
00683
00684 for (unsigned i=0; i<this->mrQuadMesh.GetNumVertices(); i++)
00685 {
00686 mPressures[i] = this->mCurrentSolution[DIM*this->mrQuadMesh.GetNumNodes() + i];
00687 }
00688 return mPressures;
00689 }
00690
00692
00694
00695 template class IncompressibleNonlinearElasticitySolver<2>;
00696 template class IncompressibleNonlinearElasticitySolver<3>;