36 #ifndef ABSTRACTCONTINUUMMECHANICSASSEMBLER_HPP_ 37 #define ABSTRACTCONTINUUMMECHANICSASSEMBLER_HPP_ 39 #include "AbstractFeAssemblerInterface.hpp" 40 #include "AbstractTetrahedralMesh.hpp" 41 #include "QuadraticMesh.hpp" 42 #include "DistributedQuadraticMesh.hpp" 43 #include "LinearBasisFunction.hpp" 44 #include "QuadraticBasisFunction.hpp" 45 #include "ReplicatableVector.hpp" 46 #include "DistributedVector.hpp" 48 #include "PetscVecTools.hpp" 49 #include "PetscMatTools.hpp" 50 #include "GaussianQuadratureRule.hpp" 86 template<
unsigned DIM,
bool CAN_ASSEMBLE_VECTOR,
bool CAN_ASSEMBLE_MATRIX>
144 c_vector<double, NUM_NODES_PER_ELEMENT>& rQuadPhi,
145 c_matrix<double, DIM, NUM_NODES_PER_ELEMENT>& rGradQuadPhi,
146 c_vector<double,DIM>& rX,
175 c_vector<double, NUM_NODES_PER_ELEMENT>& rQuadPhi,
176 c_matrix<double, DIM, NUM_NODES_PER_ELEMENT>& rGradQuadPhi,
177 c_vector<double, NUM_VERTICES_PER_ELEMENT>& rLinearPhi,
178 c_matrix<double, DIM, NUM_VERTICES_PER_ELEMENT>& rGradLinearPhi,
179 c_vector<double,DIM>& rX,
206 c_vector<double, NUM_VERTICES_PER_ELEMENT>& rLinearPhi,
207 c_matrix<double, DIM, NUM_VERTICES_PER_ELEMENT>& rGradLinearPhi,
208 c_vector<double,DIM>& rX,
238 c_vector<double, NUM_NODES_PER_ELEMENT>& rQuadPhi,
239 c_matrix<double, DIM, NUM_NODES_PER_ELEMENT>& rGradQuadPhi,
240 c_vector<double,DIM>& rX,
267 c_vector<double, NUM_VERTICES_PER_ELEMENT>& rLinearPhi,
268 c_matrix<double, DIM, NUM_VERTICES_PER_ELEMENT>& rGradLinearPhi,
269 c_vector<double,DIM>& rX,
287 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE >& rAElem,
288 c_vector<double, STENCIL_SIZE>& rBElem);
304 if ((p_quad_mesh == NULL) && (p_distributed_quad_mesh == NULL))
306 EXCEPTION(
"Continuum mechanics assemblers require a quadratic mesh");
343 template<
unsigned DIM,
bool CAN_ASSEMBLE_VECTOR,
bool CAN_ASSEMBLE_MATRIX>
351 EXCEPTION(
"Matrix to be assembled has not been set");
363 EXCEPTION(
"Vector to be assembled has not been set");
386 c_vector<double, STENCIL_SIZE> b_elem = zero_vector<double>(
STENCIL_SIZE);
415 for (
unsigned j=0; j<DIM; j++)
424 p_indices[DIM*NUM_NODES_PER_ELEMENT + i] = (DIM+1)*r_element.
GetNodeGlobalIndex(i)+DIM;
430 PetscMatTools::AddMultipleValues<STENCIL_SIZE>(this->
mMatrixToAssemble, p_indices, a_elem);
435 PetscVecTools::AddMultipleValues<STENCIL_SIZE>(this->
mVectorToAssemble, p_indices, b_elem);
441 template<
unsigned DIM,
bool CAN_ASSEMBLE_VECTOR,
bool CAN_ASSEMBLE_MATRIX>
443 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE >& rAElem,
444 c_vector<double, STENCIL_SIZE>& rBElem)
446 static c_matrix<double,DIM,DIM> jacobian;
447 static c_matrix<double,DIM,DIM> inverse_jacobian;
448 double jacobian_determinant;
463 static c_vector<double, NUM_VERTICES_PER_ELEMENT> linear_phi;
464 static c_vector<double, NUM_NODES_PER_ELEMENT> quad_phi;
465 static c_matrix<double, DIM, NUM_NODES_PER_ELEMENT> grad_quad_phi;
466 static c_matrix<double, DIM, NUM_VERTICES_PER_ELEMENT> grad_linear_phi;
468 c_vector<double,DIM> body_force;
483 c_vector<double,DIM> X = zero_vector<double>(DIM);
486 for (
unsigned j=0; j<DIM; j++)
488 X(j) += linear_phi(vertex_index)*rElement.
GetNode(vertex_index)->rGetLocation()(j);
494 c_vector<double,SPATIAL_BLOCK_SIZE_ELEMENTAL> b_spatial
496 c_vector<double,PRESSURE_BLOCK_SIZE_ELEMENTAL> b_pressure =
ComputePressureVectorTerm(linear_phi, grad_linear_phi, X, &rElement);
500 rBElem(i) += b_spatial(i)*wJ;
505 rBElem(SPATIAL_BLOCK_SIZE_ELEMENTAL + i) += b_pressure(i)*wJ;
511 c_matrix<double,SPATIAL_BLOCK_SIZE_ELEMENTAL,SPATIAL_BLOCK_SIZE_ELEMENTAL> a_spatial_spatial
514 c_matrix<double,SPATIAL_BLOCK_SIZE_ELEMENTAL,PRESSURE_BLOCK_SIZE_ELEMENTAL> a_spatial_pressure
517 c_matrix<double,PRESSURE_BLOCK_SIZE_ELEMENTAL,SPATIAL_BLOCK_SIZE_ELEMENTAL> a_pressure_spatial;
518 if (!BLOCK_SYMMETRIC_MATRIX)
524 c_matrix<double,PRESSURE_BLOCK_SIZE_ELEMENTAL,PRESSURE_BLOCK_SIZE_ELEMENTAL> a_pressure_pressure
531 rAElem(i,j) += a_spatial_spatial(i,j)*wJ;
536 rAElem(i, SPATIAL_BLOCK_SIZE_ELEMENTAL + j) += a_spatial_pressure(i,j)*wJ;
542 if (BLOCK_SYMMETRIC_MATRIX)
546 rAElem(SPATIAL_BLOCK_SIZE_ELEMENTAL + i, j) += a_spatial_pressure(j,i)*wJ;
556 rAElem(SPATIAL_BLOCK_SIZE_ELEMENTAL + i, SPATIAL_BLOCK_SIZE_ELEMENTAL + j) += a_pressure_pressure(i,j)*wJ;
563 #endif // ABSTRACTCONTINUUMMECHANICSASSEMBLER_HPP_ static const unsigned SPATIAL_BLOCK_SIZE_ELEMENTAL
bool mZeroVectorBeforeAssembly
ElementIterator GetElementIteratorBegin(bool skipDeletedElements=true)
AbstractTetrahedralMesh< DIM, DIM > * mpMesh
virtual c_matrix< double, SPATIAL_BLOCK_SIZE_ELEMENTAL, SPATIAL_BLOCK_SIZE_ELEMENTAL > ComputeSpatialSpatialMatrixTerm(c_vector< double, NUM_NODES_PER_ELEMENT > &rQuadPhi, c_matrix< double, DIM, NUM_NODES_PER_ELEMENT > &rGradQuadPhi, c_vector< double, DIM > &rX, Element< DIM, DIM > *pElement)
static const unsigned PRESSURE_BLOCK_SIZE_ELEMENTAL
virtual c_vector< double, SPATIAL_BLOCK_SIZE_ELEMENTAL > ComputeSpatialVectorTerm(c_vector< double, NUM_NODES_PER_ELEMENT > &rQuadPhi, c_matrix< double, DIM, NUM_NODES_PER_ELEMENT > &rGradQuadPhi, c_vector< double, DIM > &rX, Element< DIM, DIM > *pElement)=0
GaussianQuadratureRule< DIM > * mpQuadRule
bool mZeroMatrixBeforeAssembly
unsigned GetNodeGlobalIndex(unsigned localIndex) const
ElementIterator GetElementIteratorEnd()
virtual unsigned GetNumElements() 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)
AbstractContinuumMechanicsAssembler(AbstractTetrahedralMesh< DIM, DIM > *pMesh)
static const bool BLOCK_SYMMETRIC_MATRIX
virtual unsigned GetNumNodes() const
unsigned GetNumQuadPoints() const
virtual void GetInverseJacobianForElement(unsigned elementIndex, c_matrix< double, SPACE_DIM, ELEMENT_DIM > &rJacobian, double &rJacobianDeterminant, c_matrix< double, ELEMENT_DIM, SPACE_DIM > &rInverseJacobian) const
bool OptionExists(const std::string &rOption)
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
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 > &rReturnValue)
bool GetOwnership() const
static const unsigned STENCIL_SIZE
static void ComputeBasisFunctions(const ChastePoint< ELEMENT_DIM > &rPoint, c_vector< double,(ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2 > &rReturnValue)
virtual c_matrix< double, PRESSURE_BLOCK_SIZE_ELEMENTAL, PRESSURE_BLOCK_SIZE_ELEMENTAL > ComputePressurePressureMatrixTerm(c_vector< double, NUM_VERTICES_PER_ELEMENT > &rLinearPhi, c_matrix< double, DIM, NUM_VERTICES_PER_ELEMENT > &rGradLinearPhi, c_vector< double, DIM > &rX, Element< DIM, DIM > *pElement)
static void ComputeBasisFunctions(const ChastePoint< ELEMENT_DIM > &rPoint, c_vector< double, ELEMENT_DIM+1 > &rReturnValue)
virtual c_vector< double, PRESSURE_BLOCK_SIZE_ELEMENTAL > ComputePressureVectorTerm(c_vector< double, NUM_VERTICES_PER_ELEMENT > &rLinearPhi, c_matrix< double, DIM, NUM_VERTICES_PER_ELEMENT > &rGradLinearPhi, c_vector< double, DIM > &rX, Element< DIM, DIM > *pElement)
virtual c_matrix< double, SPATIAL_BLOCK_SIZE_ELEMENTAL, PRESSURE_BLOCK_SIZE_ELEMENTAL > ComputeSpatialPressureMatrixTerm(c_vector< double, NUM_NODES_PER_ELEMENT > &rQuadPhi, c_matrix< double, DIM, NUM_NODES_PER_ELEMENT > &rGradQuadPhi, c_vector< double, NUM_VERTICES_PER_ELEMENT > &rLinearPhi, c_matrix< double, DIM, NUM_VERTICES_PER_ELEMENT > &rGradLinearPhi, c_vector< double, DIM > &rX, Element< DIM, DIM > *pElement)
static const unsigned NUM_NODES_PER_ELEMENT
virtual ~AbstractContinuumMechanicsAssembler()
unsigned GetIndex() const
double GetWeight(unsigned index) const
static CommandLineArguments * Instance()
void AssembleOnElement(Element< DIM, DIM > &rElement, c_matrix< double, STENCIL_SIZE, STENCIL_SIZE > &rAElem, c_vector< double, STENCIL_SIZE > &rBElem)
const ChastePoint< ELEMENT_DIM > & rGetQuadPoint(unsigned index) const
static const unsigned NUM_VERTICES_PER_ELEMENT