AbstractFeObjectAssembler.hpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 #ifndef ABSTRACTFEOBJECTASSEMBLER_HPP_
00030 #define ABSTRACTFEOBJECTASSEMBLER_HPP_
00031 
00032 #include "LinearBasisFunction.hpp"
00033 #include "GaussianQuadratureRule.hpp"
00034 #include "BoundaryConditionsContainer.hpp"
00035 #include "ReplicatableVector.hpp"
00036 #include "DistributedVector.hpp"
00037 #include "HeartEventHandler.hpp"
00038 #include "AbstractTetrahedralMesh.hpp"
00039 #include "PetscTools.hpp"
00040 
00050 typedef enum InterpolationLevel_
00051 {
00052     CARDIAC = 0,
00053     NORMAL,
00054     NONLINEAR
00055 } InterpolationLevel;
00056 
00057 
00086 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL>
00087 class AbstractFeObjectAssembler : boost::noncopyable
00088 {
00089 protected:
00090 
00092     Vec mVectorToAssemble;
00093 
00095     Mat mMatrixToAssemble;
00096 
00101     bool mAssembleMatrix;
00102 
00104     bool mAssembleVector;
00105     
00107     bool mZeroVectorBeforeAssembly;
00108 
00110     bool mApplyNeummanBoundaryConditionsToVector;
00111     
00113     bool mOnlyAssembleOnSurfaceElements;
00114     
00116     PetscInt mOwnershipRangeLo;
00118     PetscInt mOwnershipRangeHi;
00119 
00121     GaussianQuadratureRule<ELEMENT_DIM>* mpQuadRule;
00122 
00124     GaussianQuadratureRule<ELEMENT_DIM-1>* mpSurfaceQuadRule;
00125 
00127     typedef LinearBasisFunction<ELEMENT_DIM> BasisFunction;
00129     typedef LinearBasisFunction<ELEMENT_DIM-1> SurfaceBasisFunction;
00130 
00131 
00136     ReplicatableVector mCurrentSolutionOrGuessReplicated;
00137 
00138 
00158     void ComputeTransformedBasisFunctionDerivatives(const ChastePoint<ELEMENT_DIM>& rPoint,
00159                                                     const c_matrix<double, ELEMENT_DIM, SPACE_DIM>& rInverseJacobian,
00160                                                     c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>& rReturnValue);
00161 
00162 
00168     void DoAssemble();
00169 
00176     void ApplyNeummanBoundaryConditionsToVector();
00177     
00183     virtual double GetCurrentSolutionOrGuessValue(unsigned nodeIndex, unsigned indexOfUnknown)
00184     {
00185         return mCurrentSolutionOrGuessReplicated[ PROBLEM_DIM*nodeIndex + indexOfUnknown];
00186     }
00187 
00188 protected:
00190     AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>* mpMesh;
00191 
00196     BoundaryConditionsContainer<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>* mpBoundaryConditions;
00197 
00203     virtual void ResetInterpolatedQuantities()
00204     {}
00205 
00214     virtual void IncrementInterpolatedQuantities(double phiI, const Node<SPACE_DIM>* pNode)
00215     {}    
00216 
00217 
00239     virtual c_matrix<double,PROBLEM_DIM*(ELEMENT_DIM+1),PROBLEM_DIM*(ELEMENT_DIM+1)> ComputeMatrixTerm(
00240         c_vector<double, ELEMENT_DIM+1>& rPhi,
00241         c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>& rGradPhi,
00242         ChastePoint<SPACE_DIM>& rX,
00243         c_vector<double,PROBLEM_DIM>& rU,
00244         c_matrix<double, PROBLEM_DIM, SPACE_DIM>& rGradU,
00245         Element<ELEMENT_DIM,SPACE_DIM>* pElement)
00246     {
00247         NEVER_REACHED;
00248         return zero_matrix<double>(PROBLEM_DIM*(ELEMENT_DIM+1),PROBLEM_DIM*(ELEMENT_DIM+1));
00249     }
00250 
00272     virtual c_vector<double,PROBLEM_DIM*(ELEMENT_DIM+1)> ComputeVectorTerm(
00273         c_vector<double, ELEMENT_DIM+1>& rPhi,
00274         c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>& rGradPhi,
00275         ChastePoint<SPACE_DIM>& rX,
00276         c_vector<double,PROBLEM_DIM>& rU,
00277         c_matrix<double, PROBLEM_DIM, SPACE_DIM>& rGradU,
00278         Element<ELEMENT_DIM,SPACE_DIM>* pElement)
00279     {
00280         NEVER_REACHED;
00281         return zero_vector<double>(PROBLEM_DIM*(ELEMENT_DIM+1));
00282     }
00283 
00301     virtual c_vector<double, PROBLEM_DIM*ELEMENT_DIM> ComputeVectorSurfaceTerm(
00302         const BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>& rSurfaceElement,
00303         c_vector<double, ELEMENT_DIM>& rPhi,
00304         ChastePoint<SPACE_DIM>& rX)
00305     {
00306         NEVER_REACHED;
00307         return zero_vector<double>(PROBLEM_DIM*ELEMENT_DIM);
00308     }
00309     
00316     virtual bool ElementAssemblyCriterion(Element<ELEMENT_DIM,SPACE_DIM>& rElement)
00317     {
00318         return true;
00319     }
00320 
00335     virtual void AssembleOnElement(Element<ELEMENT_DIM,SPACE_DIM>& rElement,
00336                                    c_matrix<double, PROBLEM_DIM*(ELEMENT_DIM+1), PROBLEM_DIM*(ELEMENT_DIM+1) >& rAElem,
00337                                    c_vector<double, PROBLEM_DIM*(ELEMENT_DIM+1)>& rBElem);
00338 
00348     virtual void AssembleOnSurfaceElement(const BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>& rSurfaceElement,
00349                                           c_vector<double, PROBLEM_DIM*ELEMENT_DIM>& rBSurfElem);
00350 
00351 
00352 public:
00359     AbstractFeObjectAssembler(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh, unsigned numQuadPoints=2);
00360 
00366     void SetApplyNeummanBoundaryConditionsToVector(BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>* pBcc);
00367 
00376     void OnlyAssembleOnSurfaceElements(bool onlyAssembleOnSurfaceElements = true)
00377     {
00378         assert(CAN_ASSEMBLE_VECTOR);
00379         mOnlyAssembleOnSurfaceElements = onlyAssembleOnSurfaceElements;
00380     }
00381     
00386     void SetMatrixToAssemble(Mat& rMatToAssemble);
00387 
00394     void SetVectorToAssemble(Vec& rVecToAssemble, bool zeroVectorBeforeAssembly);
00395 
00401     void SetCurrentSolution(Vec currentSolution);
00402 
00403 
00407     void Assemble()
00408     {
00409         mAssembleMatrix = CAN_ASSEMBLE_MATRIX;
00410         mAssembleVector = CAN_ASSEMBLE_VECTOR;
00411         DoAssemble();
00412     }
00413         
00417     void AssembleMatrix()
00418     {
00419         assert(CAN_ASSEMBLE_MATRIX);
00420         mAssembleMatrix = true;
00421         mAssembleVector = false;
00422         DoAssemble();
00423     }
00424 
00428     void AssembleVector()
00429     {
00430         assert(CAN_ASSEMBLE_VECTOR);
00431         mAssembleMatrix = false;
00432         mAssembleVector = true;
00433         DoAssemble();
00434     }
00435 
00439     virtual ~AbstractFeObjectAssembler()
00440     {
00441         delete mpQuadRule;
00442         if(mpSurfaceQuadRule)
00443         {
00444             delete mpSurfaceQuadRule;
00445         }
00446     }
00447 };
00448 
00449 
00450 
00451 
00452 
00453 
00454 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL> // class CONCRETE>
00455 AbstractFeObjectAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CAN_ASSEMBLE_VECTOR, CAN_ASSEMBLE_MATRIX, INTERPOLATION_LEVEL>::AbstractFeObjectAssembler(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh, unsigned numQuadPoints)
00456     : mVectorToAssemble(NULL),
00457       mMatrixToAssemble(NULL),
00458       mZeroVectorBeforeAssembly(true),
00459       mApplyNeummanBoundaryConditionsToVector(false),
00460       mOnlyAssembleOnSurfaceElements(false),
00461       mpMesh(pMesh),
00462       mpBoundaryConditions(NULL)
00463 {
00464     assert(pMesh);
00465     assert(numQuadPoints > 0);
00466     assert(CAN_ASSEMBLE_VECTOR || CAN_ASSEMBLE_MATRIX);
00467 
00468     mpQuadRule = new GaussianQuadratureRule<ELEMENT_DIM>(numQuadPoints);
00469 
00470     if(CAN_ASSEMBLE_VECTOR)
00471     {
00472         mpSurfaceQuadRule = new GaussianQuadratureRule<ELEMENT_DIM-1>(numQuadPoints);
00473     }
00474     else
00475     {
00476         mpSurfaceQuadRule = NULL;
00477     }
00478 }  
00479 
00480 
00481 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL> //, class CONCRETE>
00482 void AbstractFeObjectAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CAN_ASSEMBLE_VECTOR, CAN_ASSEMBLE_MATRIX, INTERPOLATION_LEVEL>::SetMatrixToAssemble(Mat& rMatToAssemble)
00483 {
00484     assert(rMatToAssemble);
00485     MatGetOwnershipRange(rMatToAssemble, &mOwnershipRangeLo, &mOwnershipRangeHi);
00486 
00487     mMatrixToAssemble = rMatToAssemble;
00488 }
00489 
00490 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL> //, class CONCRETE>
00491 void AbstractFeObjectAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CAN_ASSEMBLE_VECTOR, CAN_ASSEMBLE_MATRIX, INTERPOLATION_LEVEL>::SetVectorToAssemble(Vec& rVecToAssemble, bool zeroVectorBeforeAssembly)
00492 {
00493     assert(rVecToAssemble);
00494     VecGetOwnershipRange(rVecToAssemble, &mOwnershipRangeLo, &mOwnershipRangeHi);
00495 
00496     mVectorToAssemble = rVecToAssemble;
00497     mZeroVectorBeforeAssembly = zeroVectorBeforeAssembly;
00498 }
00499 
00500 
00501 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL> //, class CONCRETE>
00502 void AbstractFeObjectAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CAN_ASSEMBLE_VECTOR, CAN_ASSEMBLE_MATRIX, INTERPOLATION_LEVEL>::SetCurrentSolution(Vec currentSolution)
00503 {
00504     assert(currentSolution != NULL);
00505     // Replicate the current solution and store so can be used in
00506     // AssembleOnElement
00507     
00508     HeartEventHandler::BeginEvent(HeartEventHandler::COMMUNICATION);
00509     mCurrentSolutionOrGuessReplicated.ReplicatePetscVector(currentSolution);
00510     HeartEventHandler::EndEvent(HeartEventHandler::COMMUNICATION);
00511     
00512     // the AssembleOnElement type methods will determine if a current solution or
00513     // current guess exists by looking at the size of the replicated vector, so
00514     // check the size is zero if there isn't a current solution
00515     assert( mCurrentSolutionOrGuessReplicated.GetSize()>0 );
00516 }
00517 
00518 
00519 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL> //, class CONCRETE>
00520 void AbstractFeObjectAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CAN_ASSEMBLE_VECTOR, CAN_ASSEMBLE_MATRIX, INTERPOLATION_LEVEL>::DoAssemble()
00521 {
00522     HeartEventHandler::EventType assemble_event;
00523     if (mAssembleMatrix)
00524     {
00525         assemble_event = HeartEventHandler::ASSEMBLE_SYSTEM;
00526     }
00527     else
00528     {
00529         assemble_event = HeartEventHandler::ASSEMBLE_RHS;
00530     }
00531 
00532     if(mAssembleMatrix && mMatrixToAssemble==NULL)
00533     {
00534         EXCEPTION("Matrix to be assembled has not been set");
00535     }
00536     if(mAssembleVector && mVectorToAssemble==NULL)
00537     {
00538         EXCEPTION("Vector to be assembled has not been set");
00539     }
00540 
00541     // this has to be below PrepareForAssembleSystem as in that
00542     // method the odes are solved in cardiac problems
00543     HeartEventHandler::BeginEvent(assemble_event);
00544 
00545     // Zero the matrix/vector if it is to be assembled
00546     if (mAssembleVector && mZeroVectorBeforeAssembly)
00547     {
00548         PetscVecTools::Zero(mVectorToAssemble);
00549     }
00550     if (mAssembleMatrix)
00551     {
00552         PetscMatTools::Zero(mMatrixToAssemble);
00553     }
00554 
00555     const size_t STENCIL_SIZE=PROBLEM_DIM*(ELEMENT_DIM+1);
00556     c_matrix<double, STENCIL_SIZE, STENCIL_SIZE> a_elem;
00557     c_vector<double, STENCIL_SIZE> b_elem;
00558 
00560     // loop over elements
00562     if(mAssembleMatrix || (mAssembleVector && !mOnlyAssembleOnSurfaceElements))
00563     {
00564         for (typename AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ElementIterator iter = mpMesh->GetElementIteratorBegin();
00565              iter != mpMesh->GetElementIteratorEnd();
00566              ++iter)
00567         {
00568             Element<ELEMENT_DIM, SPACE_DIM>& r_element = *iter;
00569     
00570             //Test for ownership first, since it's pointless to test the criterion on something which we might know nothing about.
00571             if ( r_element.GetOwnership() == true && ElementAssemblyCriterion(r_element)==true )
00572             {
00573                 AssembleOnElement(r_element, a_elem, b_elem);
00574     
00575                 unsigned p_indices[STENCIL_SIZE];
00576                 r_element.GetStiffnessMatrixGlobalIndices(PROBLEM_DIM, p_indices);
00577     
00578                 if (mAssembleMatrix)
00579                 {
00580                     PetscMatTools::AddMultipleValues<STENCIL_SIZE>(mMatrixToAssemble, p_indices, a_elem);
00581                 }
00582     
00583                 if (mAssembleVector)
00584                 {
00585                     PetscVecTools::AddMultipleValues<STENCIL_SIZE>(mVectorToAssemble, p_indices, b_elem);
00586                 }
00587             }
00588         }
00589     }
00590 
00592     // Apply any Neumann boundary conditions
00593     //
00594     // NB. We assume that if an element has a boundary condition on any unknown there is a boundary condition
00595     // on unknown 0. This can be so for any problem by adding zero constant conditions where required
00596     // although this is a bit inefficient. Proper solution involves changing BCC to have a map of arrays
00597     // boundary conditions rather than an array of maps.
00599 
00600 
00601     if (mApplyNeummanBoundaryConditionsToVector && mAssembleVector)
00602     {
00603         HeartEventHandler::EndEvent(assemble_event);
00604         ApplyNeummanBoundaryConditionsToVector();
00605         HeartEventHandler::BeginEvent(assemble_event);
00606     }
00607     
00608     HeartEventHandler::EndEvent(assemble_event);
00609 }
00610 
00611 
00612 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL>
00613 void AbstractFeObjectAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CAN_ASSEMBLE_VECTOR, CAN_ASSEMBLE_MATRIX, INTERPOLATION_LEVEL>::SetApplyNeummanBoundaryConditionsToVector(BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>* pBcc)
00614 {
00615     assert(CAN_ASSEMBLE_VECTOR);
00616     mApplyNeummanBoundaryConditionsToVector = true;
00617     mpBoundaryConditions = pBcc;
00618 }
00619 
00620 
00621 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL>
00622 void AbstractFeObjectAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CAN_ASSEMBLE_VECTOR, CAN_ASSEMBLE_MATRIX, INTERPOLATION_LEVEL>::ApplyNeummanBoundaryConditionsToVector()
00623 {
00624     assert(mAssembleVector);
00625     assert(mpBoundaryConditions != NULL);
00626     
00627     HeartEventHandler::BeginEvent(HeartEventHandler::NEUMANN_BCS);
00628     if (mpBoundaryConditions->AnyNonZeroNeumannConditions())
00629     {
00630         typename BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::NeumannMapIterator
00631             neumann_iterator = mpBoundaryConditions->BeginNeumann();
00632         c_vector<double, PROBLEM_DIM*ELEMENT_DIM> b_surf_elem;
00633 
00634         // Iterate over defined conditions
00635         while (neumann_iterator != mpBoundaryConditions->EndNeumann())
00636         {
00637             const BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>& r_surf_element = *(neumann_iterator->first);
00638             AssembleOnSurfaceElement(r_surf_element, b_surf_elem);
00639 
00640             const size_t STENCIL_SIZE=PROBLEM_DIM*ELEMENT_DIM; // problem_dim*num_nodes_on_surface_element
00641             unsigned p_indices[STENCIL_SIZE];
00642             r_surf_element.GetStiffnessMatrixGlobalIndices(PROBLEM_DIM, p_indices);
00643             PetscVecTools::AddMultipleValues<STENCIL_SIZE>(mVectorToAssemble, p_indices, b_surf_elem);
00644             ++neumann_iterator;
00645         }
00646     }
00647     HeartEventHandler::EndEvent(HeartEventHandler::NEUMANN_BCS);
00648 }
00649 
00650 
00651 
00652 
00653 
00655 // Implementation - AssembleOnElement and smaller
00657 
00658 
00659 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL>
00660 void AbstractFeObjectAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CAN_ASSEMBLE_VECTOR, CAN_ASSEMBLE_MATRIX, INTERPOLATION_LEVEL>::ComputeTransformedBasisFunctionDerivatives(
00661         const ChastePoint<ELEMENT_DIM>& rPoint,
00662         const c_matrix<double, ELEMENT_DIM, SPACE_DIM>& rInverseJacobian,
00663         c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>& rReturnValue)
00664 {
00665     assert(ELEMENT_DIM < 4 && ELEMENT_DIM > 0);
00666     static c_matrix<double, ELEMENT_DIM, ELEMENT_DIM+1> grad_phi;
00667 
00668     LinearBasisFunction<ELEMENT_DIM>::ComputeBasisFunctionDerivatives(rPoint, grad_phi);
00669     rReturnValue = prod(trans(rInverseJacobian), grad_phi);
00670 }
00671 
00672 
00673 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL>
00674 void AbstractFeObjectAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CAN_ASSEMBLE_VECTOR, CAN_ASSEMBLE_MATRIX, INTERPOLATION_LEVEL>::AssembleOnElement(
00675     Element<ELEMENT_DIM,SPACE_DIM>& rElement,
00676     c_matrix<double, PROBLEM_DIM*(ELEMENT_DIM+1), PROBLEM_DIM*(ELEMENT_DIM+1) >& rAElem,
00677     c_vector<double, PROBLEM_DIM*(ELEMENT_DIM+1)>& rBElem)
00678 {
00684     c_matrix<double, SPACE_DIM, ELEMENT_DIM> jacobian;
00685     c_matrix<double, ELEMENT_DIM, SPACE_DIM> inverse_jacobian;
00686     double jacobian_determinant;
00687 
00688     mpMesh->GetInverseJacobianForElement(rElement.GetIndex(), jacobian, jacobian_determinant, inverse_jacobian);
00689 
00690     // With the new signature of GetInverseJacobianForElement, inverse and jacobian are returned at the same time
00691     //        // Initialise element contributions to zero
00692     //        if ( mAssembleMatrix || this->ProblemIsNonlinear() ) // don't need to construct grad_phi or grad_u in that case
00693     //        {
00694     //            this->mpMesh->GetInverseJacobianForElement(rElement.GetIndex(), inverse_jacobian);
00695     //        }
00696 
00697     if (mAssembleMatrix)
00698     {
00699         rAElem.clear();
00700     }
00701 
00702     if (mAssembleVector)
00703     {
00704         rBElem.clear();
00705     }
00706 
00707     const unsigned num_nodes = rElement.GetNumNodes();
00708 
00709     // allocate memory for the basis functions values and derivative values
00710     c_vector<double, ELEMENT_DIM+1> phi;
00711     c_matrix<double, SPACE_DIM, ELEMENT_DIM+1> grad_phi;
00712 
00713     // loop over Gauss points
00714     for (unsigned quad_index=0; quad_index < mpQuadRule->GetNumQuadPoints(); quad_index++)
00715     {
00716         const ChastePoint<ELEMENT_DIM>& quad_point = mpQuadRule->rGetQuadPoint(quad_index);
00717 
00718         BasisFunction::ComputeBasisFunctions(quad_point, phi);
00719 
00720         if ( mAssembleMatrix || INTERPOLATION_LEVEL==NONLINEAR )
00721         {
00722             ComputeTransformedBasisFunctionDerivatives(quad_point, inverse_jacobian, grad_phi);
00723         }
00724 
00725         // Location of the gauss point in the original element will be stored in x
00726         // Where applicable, u will be set to the value of the current solution at x
00727         ChastePoint<SPACE_DIM> x(0,0,0);
00728 
00729         c_vector<double,PROBLEM_DIM> u = zero_vector<double>(PROBLEM_DIM);
00730         c_matrix<double,PROBLEM_DIM,SPACE_DIM> grad_u = zero_matrix<double>(PROBLEM_DIM,SPACE_DIM);
00731 
00732         // allow the concrete version of the assembler to interpolate any
00733         // desired quantities
00734         ResetInterpolatedQuantities();
00735 
00737         // interpolation
00739         for (unsigned i=0; i<num_nodes; i++)
00740         {
00741             const Node<SPACE_DIM>* p_node = rElement.GetNode(i);
00742 
00743             if (INTERPOLATION_LEVEL != CARDIAC) // don't interpolate X if cardiac problem
00744             {
00745                 const c_vector<double, SPACE_DIM>& r_node_loc = p_node->rGetLocation();
00746                 // interpolate x
00747                 x.rGetLocation() += phi(i)*r_node_loc;
00748             }
00749 
00750             // interpolate u and grad u if a current solution or guess exists
00751             unsigned node_global_index = rElement.GetNodeGlobalIndex(i);
00752             if (mCurrentSolutionOrGuessReplicated.GetSize()>0)
00753             {
00754                 for (unsigned index_of_unknown=0; index_of_unknown<(INTERPOLATION_LEVEL!=CARDIAC ? PROBLEM_DIM : 1); index_of_unknown++)
00755                 {
00756                     // If we have a current solution (e.g. this is a dynamic problem)
00757                     // interpolate the value at this quad point
00758 
00759                     // NOTE - following assumes that, if say there are two unknowns u and v, they
00760                     // are stored in the current solution vector as
00761                     // [U1 V1 U2 V2 ... U_n V_n]
00762                     double u_at_node=GetCurrentSolutionOrGuessValue(node_global_index, index_of_unknown);
00763                     u(index_of_unknown) += phi(i)*u_at_node;
00764 
00765                     if (INTERPOLATION_LEVEL==NONLINEAR) // don't need to construct grad_phi or grad_u in other cases
00766                     {
00767                         for (unsigned j=0; j<SPACE_DIM; j++)
00768                         {
00769                             grad_u(index_of_unknown,j) += grad_phi(j,i)*u_at_node;
00770                         }
00771                     }
00772                 }
00773             }
00774 
00775             // allow the concrete version of the assembler to interpolate any
00776             // desired quantities
00777             IncrementInterpolatedQuantities(phi(i), p_node);
00778         }
00779 
00780         double wJ = jacobian_determinant * mpQuadRule->GetWeight(quad_index);
00781 
00783         // create rAElem and rBElem
00785         if (mAssembleMatrix)
00786         {
00787             noalias(rAElem) += ComputeMatrixTerm(phi, grad_phi, x, u, grad_u, &rElement) * wJ;
00788         }
00789 
00790         if (mAssembleVector)
00791         {
00792             noalias(rBElem) += ComputeVectorTerm(phi, grad_phi, x, u, grad_u, &rElement) * wJ;
00793         }
00794     }
00795 }
00796 
00797 
00798 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL>
00799 void AbstractFeObjectAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CAN_ASSEMBLE_VECTOR, CAN_ASSEMBLE_MATRIX, INTERPOLATION_LEVEL>::AssembleOnSurfaceElement(
00800             const BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>& rSurfaceElement,
00801             c_vector<double, PROBLEM_DIM*ELEMENT_DIM>& rBSurfElem)
00802 {
00803     c_vector<double, SPACE_DIM> weighted_direction;
00804     double jacobian_determinant;
00805     mpMesh->GetWeightedDirectionForBoundaryElement(rSurfaceElement.GetIndex(), weighted_direction, jacobian_determinant);
00806 
00807     rBSurfElem.clear();
00808 
00809     // allocate memory for the basis function values
00810     c_vector<double, ELEMENT_DIM>  phi;
00811 
00812     // loop over Gauss points
00813     for (unsigned quad_index=0; quad_index<mpSurfaceQuadRule->GetNumQuadPoints(); quad_index++)
00814     {
00815         const ChastePoint<ELEMENT_DIM-1>& quad_point = mpSurfaceQuadRule->rGetQuadPoint(quad_index);
00816 
00817         SurfaceBasisFunction::ComputeBasisFunctions(quad_point, phi);
00818 
00820         // interpolation
00822 
00823         // Location of the gauss point in the original element will be
00824         // stored in x
00825         ChastePoint<SPACE_DIM> x(0,0,0);
00826 
00827         this->ResetInterpolatedQuantities();
00828         for (unsigned i=0; i<rSurfaceElement.GetNumNodes(); i++)
00829         {
00830             const c_vector<double, SPACE_DIM> node_loc = rSurfaceElement.GetNode(i)->rGetLocation();
00831             x.rGetLocation() += phi(i)*node_loc;
00832 
00833             // allow the concrete version of the assembler to interpolate any
00834             // desired quantities
00835             IncrementInterpolatedQuantities(phi(i), rSurfaceElement.GetNode(i));
00836         }
00837 
00838         double wJ = jacobian_determinant * mpSurfaceQuadRule->GetWeight(quad_index);
00839 
00841         // create rAElem and rBElem
00844         noalias(rBSurfElem) += ComputeVectorSurfaceTerm(rSurfaceElement, phi, x) * wJ;
00845     }
00846 }
00847 
00848 #endif /*ABSTRACTFEOBJECTASSEMBLER_HPP_*/

Generated on Mon Apr 18 11:35:36 2011 for Chaste by  doxygen 1.5.5