AbstractContinuumMechanicsAssembler.hpp

00001 /*
00002 
00003 Copyright (c) 2005-2015, University of Oxford.
00004 All rights reserved.
00005 
00006 University of Oxford means the Chancellor, Masters and Scholars of the
00007 University of Oxford, having an administrative office at Wellington
00008 Square, Oxford OX1 2JD, UK.
00009 
00010 This file is part of Chaste.
00011 
00012 Redistribution and use in source and binary forms, with or without
00013 modification, are permitted provided that the following conditions are met:
00014  * Redistributions of source code must retain the above copyright notice,
00015    this list of conditions and the following disclaimer.
00016  * Redistributions in binary form must reproduce the above copyright notice,
00017    this list of conditions and the following disclaimer in the documentation
00018    and/or other materials provided with the distribution.
00019  * Neither the name of the University of Oxford nor the names of its
00020    contributors may be used to endorse or promote products derived from this
00021    software without specific prior written permission.
00022 
00023 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
00024 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00025 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
00026 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
00027 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
00028 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
00029 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
00030 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
00031 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
00032 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
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; //generalise to non-block symmetric matrices later (when needed maybe)
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; // assuming quadratic
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         //Check that the mesh is Quadratic
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         // In general the Jacobian for a mechanics problem is non-polynomial.
00315         // We therefore use the highest order integration rule available
00316         mpQuadRule = new GaussianQuadratureRule<DIM>(3);
00317     }
00318 
00319 
00320 //    void SetCurrentSolution(Vec currentSolution);
00321 
00322 
00326     virtual ~AbstractContinuumMechanicsAssembler()
00327     {
00328         delete mpQuadRule;
00329     }
00330 };
00331 
00332 
00334 //template<unsigned DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX>
00335 //void AbstractContinuumMechanicsAssembler<DIM,CAN_ASSEMBLE_VECTOR,CAN_ASSEMBLE_MATRIX>::SetCurrentSolution(Vec currentSolution)
00336 //{
00337 //    assert(currentSolution != NULL);
00338 //
00339 //    // Replicate the current solution and store so can be used in AssembleOnElement
00340 //    HeartEventHandler::BeginEvent(HeartEventHandler::COMMUNICATION);
00341 //    mCurrentSolutionOrGuessReplicated.ReplicatePetscVector(currentSolution);
00342 //    HeartEventHandler::EndEvent(HeartEventHandler::COMMUNICATION);
00343 //
00344 //    // The AssembleOnElement type methods will determine if a current solution or
00345 //    // current guess exists by looking at the size of the replicated vector, so
00346 //    // check the size is zero if there isn't a current solution.
00347 //    assert(mCurrentSolutionOrGuessReplicated.GetSize() > 0);
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     // Zero the matrix/vector if it is to be assembled
00379     if (this->mAssembleVector && this->mZeroVectorBeforeAssembly)
00380     {
00381         // Note PetscVecTools::Finalise(this->mVectorToAssemble); on an unused matrix
00382         // would "compress" data and make any pre-allocated entries redundant.
00383         PetscVecTools::Zero(this->mVectorToAssemble);
00384     }
00385     if (this->mAssembleMatrix && this->mZeroMatrixBeforeAssembly)
00386     {
00387         // Note PetscMatTools::Finalise(this->mMatrixToAssemble); on an unused matrix
00388         // would "compress" data and make any pre-allocated entries redundant.
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     // Loop over elements
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         // Test for ownership first, since it's pointless to test the criterion on something which we might know nothing about.
00404         if ( r_element.GetOwnership() == true  /*&& ElementAssemblyCriterion(r_element)==true*/ )
00405         {
00406             #define COVERAGE_IGNORE
00407             // note: if assemble matrix only
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             // Note that a different ordering is used for the elemental matrix compared to the global matrix.
00418             // See comments about ordering above.
00419             unsigned p_indices[STENCIL_SIZE];
00420             // Work out the mapping for spatial terms
00421             for (unsigned i=0; i<NUM_NODES_PER_ELEMENT; i++)
00422             {
00423                 for (unsigned j=0; j<DIM; j++)
00424                 {
00425                     // DIM+1 on the right-hand side here is the problem dimension
00426                     p_indices[DIM*i+j] = (DIM+1)*r_element.GetNodeGlobalIndex(i) + j;
00427                 }
00428             }
00429             // Work out the mapping for pressure terms
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     // Allocate memory for the basis functions values and derivative values
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     // Loop over Gauss points
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         // Set up basis function info
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         // interpolate X (ie physical location of this quad point).
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; // to-come: non-mixed problems
00532                 //a_pressure_spatial = ComputeSpatialPressureMatrixTerm(quad_phi, grad_quad_phi, lin_phi, grad_lin_phi, x, &rElement);
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; // to-come: non-mixed problems
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_

Generated by  doxygen 1.6.2