183 c_matrix<double, PROBLEM_DIM*NUM_CABLE_ELEMENT_NODES, PROBLEM_DIM*NUM_CABLE_ELEMENT_NODES >& rAElem,
240 if (this->mAssembleMatrix)
242 assemble_event = HeartEventHandler::ASSEMBLE_SYSTEM;
246 assemble_event = HeartEventHandler::ASSEMBLE_RHS;
249 if (this->mAssembleMatrix && this->mMatrixToAssemble==NULL)
251 EXCEPTION(
"Matrix to be assembled has not been set");
253 if (this->mAssembleVector && this->mVectorToAssemble==NULL)
255 EXCEPTION(
"Vector to be assembled has not been set");
262 if (this->mAssembleVector && this->mZeroVectorBeforeAssembly)
266 if (this->mAssembleMatrix && this->mZeroMatrixBeforeAssembly)
271 const size_t STENCIL_SIZE=PROBLEM_DIM*NUM_CABLE_ELEMENT_NODES;
272 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE> a_elem;
273 c_vector<double, STENCIL_SIZE> b_elem;
276 if (this->mAssembleMatrix || this->mAssembleVector)
279 iter != mpMesh->GetCableElementIteratorEnd();
285 if (r_element.
GetOwnership() ==
true && ElementAssemblyCriterion(r_element)==
true)
287 AssembleOnCableElement(r_element, a_elem, b_elem);
289 unsigned p_indices[STENCIL_SIZE];
292 if (this->mAssembleMatrix)
294 PetscMatTools::AddMultipleValues<STENCIL_SIZE>(this->mMatrixToAssemble, p_indices, a_elem);
297 if (this->mAssembleVector)
299 PetscVecTools::AddMultipleValues<STENCIL_SIZE>(this->mVectorToAssemble, p_indices, b_elem);
327 c_matrix<double, PROBLEM_DIM*NUM_CABLE_ELEMENT_NODES, PROBLEM_DIM*NUM_CABLE_ELEMENT_NODES >& rAElem,
328 c_vector<double, PROBLEM_DIM*NUM_CABLE_ELEMENT_NODES>& rBElem)
335 c_matrix<double, SPACE_DIM, CABLE_ELEMENT_DIM> jacobian;
336 c_matrix<double, CABLE_ELEMENT_DIM, SPACE_DIM> inverse_jacobian;
337 double jacobian_determinant;
345 if (this->mAssembleMatrix)
350 if (this->mAssembleVector)
358 c_vector<double, NUM_CABLE_ELEMENT_NODES> phi;
359 c_matrix<double, SPACE_DIM, NUM_CABLE_ELEMENT_NODES> grad_phi;
362 for (
unsigned quad_index=0; quad_index < mpCableQuadRule->GetNumQuadPoints(); quad_index++)
366 CableBasisFunction::ComputeBasisFunctions(quad_point, phi);
368 if (this->mAssembleMatrix || INTERPOLATION_LEVEL==NONLINEAR)
370 ComputeTransformedBasisFunctionDerivatives(quad_point, inverse_jacobian, grad_phi);
377 c_vector<double,PROBLEM_DIM> u = zero_vector<double>(PROBLEM_DIM);
378 c_matrix<double,PROBLEM_DIM,SPACE_DIM> grad_u = zero_matrix<double>(PROBLEM_DIM,SPACE_DIM);
381 this->ResetInterpolatedQuantities();
384 for (
unsigned i=0; i<num_nodes; i++)
388 if (INTERPOLATION_LEVEL != CARDIAC)
390 const c_vector<double, SPACE_DIM>& r_node_loc = p_node->
rGetLocation();
392 x.rGetLocation() += phi(i)*r_node_loc;
397 if (this->mCurrentSolutionOrGuessReplicated.GetSize() > 0)
399 for (
unsigned index_of_unknown=0; index_of_unknown<(INTERPOLATION_LEVEL!=CARDIAC ? PROBLEM_DIM : 1); index_of_unknown++)
409 double u_at_node = this->GetCurrentSolutionOrGuessValue(node_global_index, index_of_unknown);
410 u(index_of_unknown) += phi(i)*u_at_node;
424 this->IncrementInterpolatedQuantities(phi(i), p_node);
427 double wJ = jacobian_determinant * mpCableQuadRule->GetWeight(quad_index);
430 if (this->mAssembleMatrix)
432 noalias(rAElem) += ComputeCableMatrixTerm(phi, grad_phi, x, u, grad_u, &rElement) * wJ;
435 if (this->mAssembleVector)
437 noalias(rBElem) += ComputeCableVectorTerm(phi, grad_phi, x, u, grad_u, &rElement) * wJ;