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 #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>
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>
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>
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>
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
00506
00507
00508 HeartEventHandler::BeginEvent(HeartEventHandler::COMMUNICATION);
00509 mCurrentSolutionOrGuessReplicated.ReplicatePetscVector(currentSolution);
00510 HeartEventHandler::EndEvent(HeartEventHandler::COMMUNICATION);
00511
00512
00513
00514
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>
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
00542
00543 HeartEventHandler::BeginEvent(assemble_event);
00544
00545
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
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
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
00593
00594
00595
00596
00597
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
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;
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
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
00691
00692
00693
00694
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
00710 c_vector<double, ELEMENT_DIM+1> phi;
00711 c_matrix<double, SPACE_DIM, ELEMENT_DIM+1> grad_phi;
00712
00713
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
00726
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
00733
00734 ResetInterpolatedQuantities();
00735
00737
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)
00744 {
00745 const c_vector<double, SPACE_DIM>& r_node_loc = p_node->rGetLocation();
00746
00747 x.rGetLocation() += phi(i)*r_node_loc;
00748 }
00749
00750
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
00757
00758
00759
00760
00761
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)
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
00776
00777 IncrementInterpolatedQuantities(phi(i), p_node);
00778 }
00779
00780 double wJ = jacobian_determinant * mpQuadRule->GetWeight(quad_index);
00781
00783
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
00810 c_vector<double, ELEMENT_DIM> phi;
00811
00812
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
00822
00823
00824
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
00834
00835 IncrementInterpolatedQuantities(phi(i), rSurfaceElement.GetNode(i));
00836 }
00837
00838 double wJ = jacobian_determinant * mpSurfaceQuadRule->GetWeight(quad_index);
00839
00841
00844 noalias(rBSurfElem) += ComputeVectorSurfaceTerm(rSurfaceElement, phi, x) * wJ;
00845 }
00846 }
00847
00848 #endif