576 {
578
579#ifndef NDEBUG
582 "Topological data pointer is null");
583#endif
584
586 if (!nb_dofs)
588
589 const int nb_integration_pts = data.
getN().size1();
590
591#ifndef NDEBUG
592 const int disp_rows =
topoData->objDDisplacementAtPts.size1();
593 const int disp_cols =
topoData->objDDisplacementAtPts.size2();
594 if (disp_cols != nb_integration_pts)
596 "objDDisplacementAtPts columns (%d) != nb integration points (%d)",
597 disp_cols, nb_integration_pts);
598
601 "objDDisplacementAtPts row size is %d, expected %d", disp_rows,
603 }
604#endif
605
606 auto t_w = getFTensor0IntegrationWeight();
607 const int nb_base_functions = data.
getN().size2();
609 auto t_obj_du_gamma =
610 getFTensor1FromMat<SPACE_DIM>(
topoData->objDDisplacementAtPts);
611
613
614 auto get_ftensor1 = [](
auto &
v) {
616 &
v[0], &
v[1], &
v[2]);
617 };
618
619 for (int gg = 0; gg != nb_integration_pts; ++gg) {
620 const double a = t_w * getMeasure();
621 auto t_nf = get_ftensor1(
nF);
622
623 int bb = 0;
624 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
625 t_nf(
i) +=
a * t_row_base_fun * t_obj_du_gamma(
i);
626 ++t_nf;
627 ++t_row_base_fun;
628 }
629 for (; bb != nb_base_functions; ++bb)
630 ++t_row_base_fun;
631
632 ++t_w;
633 ++t_obj_du_gamma;
634 }
635
637}
#define FTENSOR_INDEX(DIM, I)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
boost::shared_ptr< TopologicalData > topoData
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.
VectorDouble nF
local right hand side vector