Chaste Release::3.1
ContinuumMechanicsNeumannBcsAssembler.hpp
00001 /*
00002 
00003 Copyright (c) 2005-2012, 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 CONTINUUMMECHANICSNEUMANNBCSASSEMBLER_HPP_
00037 #define CONTINUUMMECHANICSNEUMANNBCSASSEMBLER_HPP_
00038 
00039 #include "AbstractFeAssemblerInterface.hpp"
00040 #include "QuadraticMesh.hpp"
00041 #include "LinearBasisFunction.hpp"
00042 #include "QuadraticBasisFunction.hpp"
00043 #include "ReplicatableVector.hpp"
00044 #include "DistributedVector.hpp"
00045 #include "PetscTools.hpp"
00046 #include "PetscVecTools.hpp"
00047 #include "PetscMatTools.hpp"
00048 #include "GaussianQuadratureRule.hpp"
00049 #include "ContinuumMechanicsProblemDefinition.hpp"
00050 
00051 
00067 template<unsigned DIM>
00068 class ContinuumMechanicsNeumannBcsAssembler : public AbstractFeAssemblerInterface<true,false>
00069 {
00071     static const unsigned NUM_VERTICES_PER_ELEMENT = DIM;
00072 
00074     static const unsigned NUM_NODES_PER_ELEMENT = DIM*(DIM+1)/2; // assuming quadratic
00075 
00077     static const unsigned SPATIAL_BLOCK_SIZE_ELEMENTAL = DIM*NUM_NODES_PER_ELEMENT;
00079     static const unsigned PRESSURE_BLOCK_SIZE_ELEMENTAL = NUM_VERTICES_PER_ELEMENT;
00080 
00082     static const unsigned STENCIL_SIZE = DIM*NUM_NODES_PER_ELEMENT + NUM_VERTICES_PER_ELEMENT;
00083 
00084 protected:
00086     QuadraticMesh<DIM>* mpMesh;
00087 
00089     ContinuumMechanicsProblemDefinition<DIM>* mpProblemDefinition;
00090 
00092     GaussianQuadratureRule<DIM-1>* mpQuadRule;
00093 
00098     void DoAssemble();
00099 
00100 
00111     void AssembleOnBoundaryElement(BoundaryElement<DIM-1,DIM>& rElement,
00112                                    c_vector<double, STENCIL_SIZE>& rBElem,
00113                                    unsigned boundaryConditionIndex);
00114 
00115 public:
00120     ContinuumMechanicsNeumannBcsAssembler(QuadraticMesh<DIM>* pMesh,
00121                                           ContinuumMechanicsProblemDefinition<DIM>* pProblemDefinition,
00122                                           unsigned numQuadPoints = 2)
00123         : AbstractFeAssemblerInterface<true,false>(),
00124           mpMesh(pMesh),
00125           mpProblemDefinition(pProblemDefinition)
00126     {
00127         assert(pMesh);
00128         assert(pProblemDefinition);
00129         mpQuadRule = new GaussianQuadratureRule<DIM-1>(numQuadPoints);
00130     }
00131 
00132 
00136     virtual ~ContinuumMechanicsNeumannBcsAssembler()
00137     {
00138         delete mpQuadRule;
00139     }
00140 };
00141 
00142 
00143 template<unsigned DIM>
00144 void ContinuumMechanicsNeumannBcsAssembler<DIM>::DoAssemble()
00145 {
00146     if(this->mVectorToAssemble==NULL)
00147     {
00148         EXCEPTION("Vector to be assembled has not been set");
00149     }
00150 
00151     if( PetscVecTools::GetSize(this->mVectorToAssemble) != (DIM+1)*mpMesh->GetNumNodes() )
00152     {
00153         EXCEPTION("Vector provided to be assembled has size " << PetscVecTools::GetSize(this->mVectorToAssemble) << ", not expected size of " << (DIM+1)*mpMesh->GetNumNodes() << " ((dim+1)*num_nodes)");
00154     }
00155 
00156     // Zero the matrix/vector if it is to be assembled
00157     if (this->mZeroVectorBeforeAssembly)
00158     {
00159         PetscVecTools::Zero(this->mVectorToAssemble);
00160     }
00161 
00162 
00163     if (mpProblemDefinition->GetTractionBoundaryConditionType() != NO_TRACTIONS)
00164     {
00165         c_vector<double, STENCIL_SIZE> b_elem = zero_vector<double>(STENCIL_SIZE);
00166 
00167         for (unsigned bc_index=0; bc_index<mpProblemDefinition->rGetTractionBoundaryElements().size(); bc_index++)
00168         {
00169             BoundaryElement<DIM-1,DIM>& r_boundary_element = *(mpProblemDefinition->rGetTractionBoundaryElements()[bc_index]);
00170             AssembleOnBoundaryElement(r_boundary_element, b_elem, bc_index);
00171 
00172             unsigned p_indices[STENCIL_SIZE];
00173             for (unsigned i=0; i<NUM_NODES_PER_ELEMENT; i++)
00174             {
00175                 for (unsigned j=0; j<DIM; j++)
00176                 {
00177                     p_indices[DIM*i+j] = (DIM+1)*r_boundary_element.GetNodeGlobalIndex(i) + j;
00178                 }
00179             }
00180             // Note: The pressure block of b_elem will be zero, but this bit still needs to be
00181             // set to avoid memory leaks.
00182             for (unsigned i=0; i<DIM /*vertices per boundary elem */; i++)
00183             {
00184                 p_indices[DIM*NUM_NODES_PER_ELEMENT + i] = (DIM+1)*r_boundary_element.GetNodeGlobalIndex(i)+DIM;
00185             }
00186 
00187             PetscVecTools::AddMultipleValues<STENCIL_SIZE>(this->mVectorToAssemble, p_indices, b_elem);
00188         }
00189     }
00190 }
00191 
00192 
00193 
00194 template<unsigned DIM>
00195 void ContinuumMechanicsNeumannBcsAssembler<DIM>::AssembleOnBoundaryElement(BoundaryElement<DIM-1,DIM>& rBoundaryElement,
00196                                                                            c_vector<double,STENCIL_SIZE>& rBelem,
00197                                                                            unsigned boundaryConditionIndex)
00198 {
00199     rBelem.clear();
00200 
00201     c_vector<double, DIM> weighted_direction;
00202     double jacobian_determinant;
00203     mpMesh->GetWeightedDirectionForBoundaryElement(rBoundaryElement.GetIndex(), weighted_direction, jacobian_determinant);
00204 
00205     c_vector<double,NUM_NODES_PER_ELEMENT> phi;
00206 
00207     for (unsigned quad_index=0; quad_index<mpQuadRule->GetNumQuadPoints(); quad_index++)
00208     {
00209         double wJ = jacobian_determinant * mpQuadRule->GetWeight(quad_index);
00210         const ChastePoint<DIM-1>& quad_point = mpQuadRule->rGetQuadPoint(quad_index);
00211         QuadraticBasisFunction<DIM-1>::ComputeBasisFunctions(quad_point, phi);
00212 
00213         c_vector<double,DIM> traction = zero_vector<double>(DIM);
00214         switch (mpProblemDefinition->GetTractionBoundaryConditionType())
00215         {
00216             case ELEMENTWISE_TRACTION:
00217             {
00218                 traction = mpProblemDefinition->rGetElementwiseTractions()[boundaryConditionIndex];
00219                 break;
00220             }
00221             default:
00222                 // Functional traction not implemented yet..
00223                 NEVER_REACHED;
00224         }
00225 
00226         for (unsigned index=0; index<NUM_NODES_PER_ELEMENT*DIM; index++)
00227         {
00228             unsigned spatial_dim = index%DIM;
00229             unsigned node_index = (index-spatial_dim)/DIM;
00230 
00231             assert(node_index < NUM_NODES_PER_ELEMENT);
00232 
00233             rBelem(index) += traction(spatial_dim) * phi(node_index) * wJ;
00234         }
00235     }
00236 }
00237 
00238 
00239 #endif // CONTINUUMMECHANICSNEUMANNBCSASSEMBLER_HPP_