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 #ifndef ABSTRACTCONTINUUMMECHANICSASSEMBLER_HPP_
00037 #define ABSTRACTCONTINUUMMECHANICSASSEMBLER_HPP_
00038
00039 #include "AbstractFeAssemblerInterface.hpp"
00040 #include "AbstractTetrahedralMesh.hpp"
00041 #include "QuadraticMesh.hpp"
00042 #include "DistributedQuadraticMesh.hpp"
00043 #include "LinearBasisFunction.hpp"
00044 #include "QuadraticBasisFunction.hpp"
00045 #include "ReplicatableVector.hpp"
00046 #include "DistributedVector.hpp"
00047 #include "PetscTools.hpp"
00048 #include "PetscVecTools.hpp"
00049 #include "PetscMatTools.hpp"
00050 #include "GaussianQuadratureRule.hpp"
00051
00052
00086 template<unsigned DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX>
00087 class AbstractContinuumMechanicsAssembler : public AbstractFeAssemblerInterface<CAN_ASSEMBLE_VECTOR,CAN_ASSEMBLE_MATRIX>
00088 {
00092 static const bool BLOCK_SYMMETRIC_MATRIX = true;
00093
00095 static const unsigned NUM_VERTICES_PER_ELEMENT = DIM+1;
00096
00098 static const unsigned NUM_NODES_PER_ELEMENT = (DIM+1)*(DIM+2)/2;
00099
00101 static const unsigned SPATIAL_BLOCK_SIZE_ELEMENTAL = DIM*NUM_NODES_PER_ELEMENT;
00103 static const unsigned PRESSURE_BLOCK_SIZE_ELEMENTAL = NUM_VERTICES_PER_ELEMENT;
00104
00106 static const unsigned STENCIL_SIZE = DIM*NUM_NODES_PER_ELEMENT + NUM_VERTICES_PER_ELEMENT;
00107
00108 protected:
00110 AbstractTetrahedralMesh<DIM,DIM>* mpMesh;
00111
00113 GaussianQuadratureRule<DIM>* mpQuadRule;
00114
00121 void DoAssemble();
00122
00123
00143 virtual c_matrix<double,SPATIAL_BLOCK_SIZE_ELEMENTAL,SPATIAL_BLOCK_SIZE_ELEMENTAL> ComputeSpatialSpatialMatrixTerm(
00144 c_vector<double, NUM_NODES_PER_ELEMENT>& rQuadPhi,
00145 c_matrix<double, DIM, NUM_NODES_PER_ELEMENT>& rGradQuadPhi,
00146 c_vector<double,DIM>& rX,
00147 Element<DIM,DIM>* pElement)
00148 {
00149 return zero_matrix<double>(SPATIAL_BLOCK_SIZE_ELEMENTAL,SPATIAL_BLOCK_SIZE_ELEMENTAL);
00150 }
00151
00174 virtual c_matrix<double,SPATIAL_BLOCK_SIZE_ELEMENTAL,PRESSURE_BLOCK_SIZE_ELEMENTAL> ComputeSpatialPressureMatrixTerm(
00175 c_vector<double, NUM_NODES_PER_ELEMENT>& rQuadPhi,
00176 c_matrix<double, DIM, NUM_NODES_PER_ELEMENT>& rGradQuadPhi,
00177 c_vector<double, NUM_VERTICES_PER_ELEMENT>& rLinearPhi,
00178 c_matrix<double, DIM, NUM_VERTICES_PER_ELEMENT>& rGradLinearPhi,
00179 c_vector<double,DIM>& rX,
00180 Element<DIM,DIM>* pElement)
00181 {
00182 return zero_matrix<double>(SPATIAL_BLOCK_SIZE_ELEMENTAL,PRESSURE_BLOCK_SIZE_ELEMENTAL);
00183 }
00184
00185
00205 virtual c_matrix<double,PRESSURE_BLOCK_SIZE_ELEMENTAL,PRESSURE_BLOCK_SIZE_ELEMENTAL> ComputePressurePressureMatrixTerm(
00206 c_vector<double, NUM_VERTICES_PER_ELEMENT>& rLinearPhi,
00207 c_matrix<double, DIM, NUM_VERTICES_PER_ELEMENT>& rGradLinearPhi,
00208 c_vector<double,DIM>& rX,
00209 Element<DIM,DIM>* pElement)
00210 {
00211 return zero_matrix<double>(PRESSURE_BLOCK_SIZE_ELEMENTAL,PRESSURE_BLOCK_SIZE_ELEMENTAL);
00212 }
00213
00214
00237 virtual c_vector<double,SPATIAL_BLOCK_SIZE_ELEMENTAL> ComputeSpatialVectorTerm(
00238 c_vector<double, NUM_NODES_PER_ELEMENT>& rQuadPhi,
00239 c_matrix<double, DIM, NUM_NODES_PER_ELEMENT>& rGradQuadPhi,
00240 c_vector<double,DIM>& rX,
00241 Element<DIM,DIM>* pElement)
00242 {
00243 return zero_vector<double>(SPATIAL_BLOCK_SIZE_ELEMENTAL);
00244 }
00245
00246
00269 virtual c_vector<double,PRESSURE_BLOCK_SIZE_ELEMENTAL> ComputePressureVectorTerm(
00270 c_vector<double, NUM_VERTICES_PER_ELEMENT>& rLinearPhi,
00271 c_matrix<double, DIM, NUM_VERTICES_PER_ELEMENT>& rGradLinearPhi,
00272 c_vector<double,DIM>& rX,
00273 Element<DIM,DIM>* pElement)
00274 {
00275 return zero_vector<double>(PRESSURE_BLOCK_SIZE_ELEMENTAL);
00276 }
00277
00278
00279
00291 void AssembleOnElement(Element<DIM, DIM>& rElement,
00292 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE >& rAElem,
00293 c_vector<double, STENCIL_SIZE>& rBElem);
00294
00295 public:
00299 AbstractContinuumMechanicsAssembler(AbstractTetrahedralMesh<DIM, DIM>* pMesh)
00300 : AbstractFeAssemblerInterface<CAN_ASSEMBLE_VECTOR,CAN_ASSEMBLE_MATRIX>(),
00301 mpMesh(pMesh)
00302 {
00303 assert(pMesh);
00304
00305
00306 QuadraticMesh<DIM>* p_quad_mesh = dynamic_cast<QuadraticMesh<DIM>* >(pMesh);
00307 DistributedQuadraticMesh<DIM>* p_distributed_quad_mesh = dynamic_cast<DistributedQuadraticMesh<DIM>* >(pMesh);
00308
00309 if(p_quad_mesh == NULL && p_distributed_quad_mesh == NULL)
00310 {
00311 EXCEPTION("Continuum mechanics assemblers require a quadratic mesh");
00312 }
00313
00314
00315
00316 mpQuadRule = new GaussianQuadratureRule<DIM>(3);
00317 }
00318
00319
00320
00321
00322
00326 virtual ~AbstractContinuumMechanicsAssembler()
00327 {
00328 delete mpQuadRule;
00329 }
00330 };
00331
00332
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350 template<unsigned DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX>
00351 void AbstractContinuumMechanicsAssembler<DIM,CAN_ASSEMBLE_VECTOR,CAN_ASSEMBLE_MATRIX>::DoAssemble()
00352 {
00353 assert(this->mAssembleMatrix || this->mAssembleVector);
00354 if (this->mAssembleMatrix)
00355 {
00356 if(this->mMatrixToAssemble==NULL)
00357 {
00358 EXCEPTION("Matrix to be assembled has not been set");
00359 }
00360 if( PetscMatTools::GetSize(this->mMatrixToAssemble) != (DIM+1)*mpMesh->GetNumNodes() )
00361 {
00362 EXCEPTION("Matrix provided to be assembled has size " << PetscMatTools::GetSize(this->mMatrixToAssemble) << ", not expected size of " << (DIM+1)*mpMesh->GetNumNodes() << " ((dim+1)*num_nodes)");
00363 }
00364 }
00365
00366 if (this->mAssembleVector)
00367 {
00368 if(this->mVectorToAssemble==NULL)
00369 {
00370 EXCEPTION("Vector to be assembled has not been set");
00371 }
00372 if( PetscVecTools::GetSize(this->mVectorToAssemble) != (DIM+1)*mpMesh->GetNumNodes() )
00373 {
00374 EXCEPTION("Vector provided to be assembled has size " << PetscVecTools::GetSize(this->mVectorToAssemble) << ", not expected size of " << (DIM+1)*mpMesh->GetNumNodes() << " ((dim+1)*num_nodes)");
00375 }
00376 }
00377
00378
00379 if (this->mAssembleVector && this->mZeroVectorBeforeAssembly)
00380 {
00381
00382
00383 PetscVecTools::Zero(this->mVectorToAssemble);
00384 }
00385 if (this->mAssembleMatrix && this->mZeroMatrixBeforeAssembly)
00386 {
00387
00388
00389 PetscMatTools::Zero(this->mMatrixToAssemble);
00390 }
00391
00392 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE> a_elem = zero_matrix<double>(STENCIL_SIZE,STENCIL_SIZE);
00393 c_vector<double, STENCIL_SIZE> b_elem = zero_vector<double>(STENCIL_SIZE);
00394
00395
00396
00397 for (typename AbstractTetrahedralMesh<DIM, DIM>::ElementIterator iter = mpMesh->GetElementIteratorBegin();
00398 iter != mpMesh->GetElementIteratorEnd();
00399 ++iter)
00400 {
00401 Element<DIM, DIM>& r_element = *iter;
00402
00403
00404 if ( r_element.GetOwnership() == true )
00405 {
00406 #define COVERAGE_IGNORE
00407
00408 if(CommandLineArguments::Instance()->OptionExists("-mech_very_verbose") && this->mAssembleMatrix)
00409 {
00410 std::cout << "\r[" << PetscTools::GetMyRank() << "]: Element " << r_element.GetIndex() << " of " << mpMesh->GetNumElements() << std::flush;
00411 }
00412 #undef COVERAGE_IGNORE
00413
00414
00415 AssembleOnElement(r_element, a_elem, b_elem);
00416
00417
00418
00419 unsigned p_indices[STENCIL_SIZE];
00420
00421 for (unsigned i=0; i<NUM_NODES_PER_ELEMENT; i++)
00422 {
00423 for (unsigned j=0; j<DIM; j++)
00424 {
00425
00426 p_indices[DIM*i+j] = (DIM+1)*r_element.GetNodeGlobalIndex(i) + j;
00427 }
00428 }
00429
00430 for (unsigned i=0; i<NUM_VERTICES_PER_ELEMENT; i++)
00431 {
00432 p_indices[DIM*NUM_NODES_PER_ELEMENT + i] = (DIM+1)*r_element.GetNodeGlobalIndex(i)+DIM;
00433 }
00434
00435
00436 if (this->mAssembleMatrix)
00437 {
00438 PetscMatTools::AddMultipleValues<STENCIL_SIZE>(this->mMatrixToAssemble, p_indices, a_elem);
00439 }
00440
00441 if (this->mAssembleVector)
00442 {
00443 PetscVecTools::AddMultipleValues<STENCIL_SIZE>(this->mVectorToAssemble, p_indices, b_elem);
00444 }
00445 }
00446 }
00447 }
00448
00449 template<unsigned DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX>
00450 void AbstractContinuumMechanicsAssembler<DIM,CAN_ASSEMBLE_VECTOR,CAN_ASSEMBLE_MATRIX>::AssembleOnElement(Element<DIM, DIM>& rElement,
00451 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE >& rAElem,
00452 c_vector<double, STENCIL_SIZE>& rBElem)
00453 {
00454 static c_matrix<double,DIM,DIM> jacobian;
00455 static c_matrix<double,DIM,DIM> inverse_jacobian;
00456 double jacobian_determinant;
00457
00458 mpMesh->GetInverseJacobianForElement(rElement.GetIndex(), jacobian, jacobian_determinant, inverse_jacobian);
00459
00460 if (this->mAssembleMatrix)
00461 {
00462 rAElem.clear();
00463 }
00464
00465 if (this->mAssembleVector)
00466 {
00467 rBElem.clear();
00468 }
00469
00470
00471
00472 static c_vector<double, NUM_VERTICES_PER_ELEMENT> linear_phi;
00473 static c_vector<double, NUM_NODES_PER_ELEMENT> quad_phi;
00474 static c_matrix<double, DIM, NUM_NODES_PER_ELEMENT> grad_quad_phi;
00475 static c_matrix<double, DIM, NUM_VERTICES_PER_ELEMENT> grad_linear_phi;
00476
00477 c_vector<double,DIM> body_force;
00478
00479
00480 for (unsigned quadrature_index=0; quadrature_index < mpQuadRule->GetNumQuadPoints(); quadrature_index++)
00481 {
00482 double wJ = jacobian_determinant * mpQuadRule->GetWeight(quadrature_index);
00483 const ChastePoint<DIM>& quadrature_point = mpQuadRule->rGetQuadPoint(quadrature_index);
00484
00485
00486 LinearBasisFunction<DIM>::ComputeBasisFunctions(quadrature_point, linear_phi);
00487 QuadraticBasisFunction<DIM>::ComputeBasisFunctions(quadrature_point, quad_phi);
00488 QuadraticBasisFunction<DIM>::ComputeTransformedBasisFunctionDerivatives(quadrature_point, inverse_jacobian, grad_quad_phi);
00489 LinearBasisFunction<DIM>::ComputeTransformedBasisFunctionDerivatives(quadrature_point, inverse_jacobian, grad_linear_phi);
00490
00491
00492 c_vector<double,DIM> X = zero_vector<double>(DIM);
00493 for (unsigned vertex_index=0; vertex_index<NUM_VERTICES_PER_ELEMENT; vertex_index++)
00494 {
00495 for(unsigned j=0; j<DIM; j++)
00496 {
00497 X(j) += linear_phi(vertex_index)*rElement.GetNode(vertex_index)->rGetLocation()(j);
00498 }
00499 }
00500
00501 if(this->mAssembleVector)
00502 {
00503 c_vector<double,SPATIAL_BLOCK_SIZE_ELEMENTAL> b_spatial
00504 = ComputeSpatialVectorTerm(quad_phi, grad_quad_phi, X, &rElement);
00505 c_vector<double,PRESSURE_BLOCK_SIZE_ELEMENTAL> b_pressure = ComputePressureVectorTerm(linear_phi, grad_linear_phi, X, &rElement);
00506
00507 for (unsigned i=0; i<SPATIAL_BLOCK_SIZE_ELEMENTAL; i++)
00508 {
00509 rBElem(i) += b_spatial(i)*wJ;
00510 }
00511
00512
00513 for (unsigned i=0; i<PRESSURE_BLOCK_SIZE_ELEMENTAL; i++)
00514 {
00515 rBElem(SPATIAL_BLOCK_SIZE_ELEMENTAL + i) += b_pressure(i)*wJ;
00516 }
00517 }
00518
00519
00520 if(this->mAssembleMatrix)
00521 {
00522 c_matrix<double,SPATIAL_BLOCK_SIZE_ELEMENTAL,SPATIAL_BLOCK_SIZE_ELEMENTAL> a_spatial_spatial
00523 = ComputeSpatialSpatialMatrixTerm(quad_phi, grad_quad_phi, X, &rElement);
00524
00525 c_matrix<double,SPATIAL_BLOCK_SIZE_ELEMENTAL,PRESSURE_BLOCK_SIZE_ELEMENTAL> a_spatial_pressure
00526 = ComputeSpatialPressureMatrixTerm(quad_phi, grad_quad_phi, linear_phi, grad_linear_phi, X, &rElement);
00527
00528 c_matrix<double,PRESSURE_BLOCK_SIZE_ELEMENTAL,SPATIAL_BLOCK_SIZE_ELEMENTAL> a_pressure_spatial;
00529 if(!BLOCK_SYMMETRIC_MATRIX)
00530 {
00531 NEVER_REACHED;
00532
00533 }
00534
00535 c_matrix<double,PRESSURE_BLOCK_SIZE_ELEMENTAL,PRESSURE_BLOCK_SIZE_ELEMENTAL> a_pressure_pressure
00536 = ComputePressurePressureMatrixTerm(linear_phi, grad_linear_phi, X, &rElement);
00537
00538 for (unsigned i=0; i<SPATIAL_BLOCK_SIZE_ELEMENTAL; i++)
00539 {
00540 for(unsigned j=0; j<SPATIAL_BLOCK_SIZE_ELEMENTAL; j++)
00541 {
00542 rAElem(i,j) += a_spatial_spatial(i,j)*wJ;
00543 }
00544
00545 for(unsigned j=0; j<PRESSURE_BLOCK_SIZE_ELEMENTAL; j++)
00546 {
00547 rAElem(i, SPATIAL_BLOCK_SIZE_ELEMENTAL + j) += a_spatial_pressure(i,j)*wJ;
00548 }
00549 }
00550
00551 for(unsigned i=0; i<PRESSURE_BLOCK_SIZE_ELEMENTAL; i++)
00552 {
00553 if(BLOCK_SYMMETRIC_MATRIX)
00554 {
00555 for(unsigned j=0; j<SPATIAL_BLOCK_SIZE_ELEMENTAL; j++)
00556 {
00557 rAElem(SPATIAL_BLOCK_SIZE_ELEMENTAL + i, j) += a_spatial_pressure(j,i)*wJ;
00558 }
00559 }
00560 else
00561 {
00562 NEVER_REACHED;
00563 }
00564
00565 for(unsigned j=0; j<PRESSURE_BLOCK_SIZE_ELEMENTAL; j++)
00566 {
00567 rAElem(SPATIAL_BLOCK_SIZE_ELEMENTAL + i, SPATIAL_BLOCK_SIZE_ELEMENTAL + j) += a_pressure_pressure(i,j)*wJ;
00568 }
00569 }
00570 }
00571 }
00572 }
00573
00574
00575 #endif // ABSTRACTCONTINUUMMECHANICSASSEMBLER_HPP_