Chaste Release::3.1
|
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_