157 if (this->mVectorToAssemble==NULL)
159 EXCEPTION(
"Vector to be assembled has not been set");
164 EXCEPTION(
"Vector provided to be assembled has size " <<
PetscVecTools::GetSize(this->mVectorToAssemble) <<
", not expected size of " << (DIM+1)*mpMesh->GetNumNodes() <<
" ((dim+1)*num_nodes)");
168 if (this->mZeroVectorBeforeAssembly)
174 if (mpProblemDefinition->GetTractionBoundaryConditionType() != NO_TRACTIONS)
176 c_vector<double, STENCIL_SIZE> b_elem = zero_vector<double>(STENCIL_SIZE);
178 for (
unsigned bc_index=0; bc_index<mpProblemDefinition->rGetTractionBoundaryElements().size(); bc_index++)
180 BoundaryElement<DIM-1,DIM>& r_boundary_element = *(mpProblemDefinition->rGetTractionBoundaryElements()[bc_index]);
181 AssembleOnBoundaryElement(r_boundary_element, b_elem, bc_index);
183 unsigned p_indices[STENCIL_SIZE];
184 for (
unsigned i=0; i<NUM_NODES_PER_ELEMENT; i++)
186 for (
unsigned j=0; j<DIM; j++)
188 p_indices[DIM*i+j] = (DIM+1)*r_boundary_element.GetNodeGlobalIndex(i) + j;
193 for (
unsigned i=0; i<DIM ; i++)
195 p_indices[DIM*NUM_NODES_PER_ELEMENT + i] = (DIM+1)*r_boundary_element.GetNodeGlobalIndex(i)+DIM;
198 PetscVecTools::AddMultipleValues<STENCIL_SIZE>(this->mVectorToAssemble, p_indices, b_elem);
205 c_vector<double,STENCIL_SIZE>& rBelem,
206 unsigned boundaryConditionIndex)
210 c_vector<double, DIM> weighted_direction;
211 double jacobian_determinant;
212 mpMesh->GetWeightedDirectionForBoundaryElement(rBoundaryElement.
GetIndex(), weighted_direction, jacobian_determinant);
214 c_vector<double,NUM_NODES_PER_ELEMENT> phi;
216 for (
unsigned quad_index=0; quad_index<mpQuadRule->GetNumQuadPoints(); quad_index++)
218 double wJ = jacobian_determinant * mpQuadRule->GetWeight(quad_index);
219 const ChastePoint<DIM-1>& quad_point = mpQuadRule->rGetQuadPoint(quad_index);
222 c_vector<double,DIM> traction = zero_vector<double>(DIM);
223 switch (mpProblemDefinition->GetTractionBoundaryConditionType())
225 case ELEMENTWISE_TRACTION:
227 traction = mpProblemDefinition->rGetElementwiseTractions()[boundaryConditionIndex];
235 for (
unsigned index=0; index<NUM_NODES_PER_ELEMENT*DIM; index++)
237 unsigned spatial_dim = index%DIM;
238 unsigned node_index = (index-spatial_dim)/DIM;
240 assert(node_index < NUM_NODES_PER_ELEMENT);
242 rBelem(index) += traction(spatial_dim) * phi(node_index) * wJ;