26 doEntities{true, true, true, true, true, true,
27 true, true, true, true, true, true},
29 doVertices(doEntities[MBVERTEX]), doEdges(doEntities[MBEDGE]),
30 doQuads(doEntities[MBQUAD]), doTris(doEntities[MBTRI]),
31 doTets(doEntities[MBTET]), doPrisms(doEntities[MBPRISM]) {
46 [&](boost::ptr_vector<EntitiesFieldData::EntData> &row_ent_data,
50 for (
EntityType col_type = low_type; col_type != hi_type; ++col_type) {
52 for (
size_t SS = 0; SS != col_ent_data.size(); SS++) {
53 if (col_ent_data[SS].getFieldData().size())
54 CHKERR doWork(ss, SS, row_type, col_type, row_ent_data[ss],
61 auto do_row_entity = [&](
const EntityType type) {
64 for (
size_t ss = 0; ss != row_ent_data.size(); ++ss) {
66 CHKERR do_col_entity(row_ent_data, ss, type, MBVERTEX, type);
74 CHKERR do_col_entity(row_ent_data, ss, type,
80 for (
EntityType row_type = MBVERTEX; row_type != MBMAXTYPE; ++row_type) {
82 CHKERR do_row_entity(row_type);
93 return opLhs<true>(row_data, col_data);
95 return opLhs<false>(row_data, col_data);
98template <
bool ErrorIfNoBase>
101 const std::array<bool, MBMAXTYPE> &do_entities) {
104 auto do_entity = [&](
auto type) {
108 const size_t size = ent_data.size();
109 for (
size_t ss = 0; ss != size; ++ss) {
111 auto &side_data = ent_data[ss];
114 if (side_data.getFieldData().size() &&
115 (side_data.getBase() ==
NOBASE ||
117 for (VectorDofs::iterator it = side_data.getFieldDofs().begin();
118 it != side_data.getFieldDofs().end(); it++)
119 if ((*it) && (*it)->getActive())
130 for (
EntityType row_type = MBVERTEX; row_type != MBMAXTYPE; ++row_type) {
131 if (do_entities[row_type]) {
132 CHKERR do_entity(row_type);
141 const bool error_if_no_base) {
142 if (error_if_no_base)
153 auto A = getFTensor2FromMat<3, 3>(jac_data);
154 int nb_gauss_pts = jac_data.size2();
155 det_data.resize(nb_gauss_pts,
false);
156 inv_jac_data.resize(3, nb_gauss_pts,
false);
158 auto I = getFTensor2FromMat<3, 3>(inv_jac_data);
159 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
176 if (diff_n.data().size()) {
177 const int nb_base_functions = diff_n.size2() / 3;
178 const int nb_gauss_pts = diff_n.size1();
179 diffNinvJac.resize(diff_n.size1(), diff_n.size2(),
false);
181 double *t_diff_n_ptr = &*diff_n.data().begin();
183 t_diff_n_ptr, &t_diff_n_ptr[1], &t_diff_n_ptr[2]);
184 double *t_inv_n_ptr = &*
diffNinvJac.data().begin();
186 t_inv_n_ptr, &t_inv_n_ptr[1], &t_inv_n_ptr[2]);
188 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
189 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
210 CHKERR transform_base(*(
m.second));
215 CHKERR transform_base(*ptr);
226 if (type == MBVERTEX)
233 const unsigned int nb_gauss_pts = data.
getDiffN(base).size1();
234 const unsigned int nb_base_functions = data.
getDiffN(base).size2() / 9;
235 if (!nb_base_functions)
243 t_inv_diff_n_ptr, &t_inv_diff_n_ptr[
HVEC0_1],
252 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
253 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
270 if (CN::Dimension(type) > 1) {
276 const unsigned int nb_base_functions = data.
getN(base).size2() / 3;
277 if (!nb_base_functions)
280 const unsigned int nb_gauss_pts = data.
getN(base).size1();
283 piolaN.resize(nb_gauss_pts, data.
getN(base).size2(),
false);
284 if (data.
getN(base).size2() > 0) {
286 double *t_transformed_n_ptr = &*
piolaN.data().begin();
289 &t_transformed_n_ptr[
HVEC1], &t_transformed_n_ptr[
HVEC2]);
290 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
291 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
292 t_transformed_n(
i) =
a * (
tJac(
i,
k) * t_n(
k));
301 if (data.
getDiffN(base).size2() > 0) {
303 double *t_transformed_diff_n_ptr = &*
piolaDiffN.data().begin();
305 t_transformed_diff_n(t_transformed_diff_n_ptr,
306 &t_transformed_diff_n_ptr[
HVEC0_1],
307 &t_transformed_diff_n_ptr[
HVEC0_2],
308 &t_transformed_diff_n_ptr[
HVEC1_0],
309 &t_transformed_diff_n_ptr[
HVEC1_1],
310 &t_transformed_diff_n_ptr[
HVEC1_2],
311 &t_transformed_diff_n_ptr[
HVEC2_0],
312 &t_transformed_diff_n_ptr[
HVEC2_1],
313 &t_transformed_diff_n_ptr[
HVEC2_2]);
314 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
315 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
316 t_transformed_diff_n(
i,
k) =
a *
tJac(
i,
j) * t_diff_n(
j,
k);
318 ++t_transformed_diff_n;
334 if (type == MBVERTEX)
341 const unsigned int nb_base_functions = data.
getN(base).size2() / 3;
342 if (!nb_base_functions)
345 const unsigned int nb_gauss_pts = data.
getN(base).size1();
346 piolaN.resize(nb_gauss_pts, data.
getN(base).size2(),
false);
350 double *t_transformed_n_ptr = &*
piolaN.data().begin();
352 t_transformed_n_ptr, &t_transformed_n_ptr[
HVEC1],
353 &t_transformed_n_ptr[
HVEC2]);
355 double *t_transformed_diff_n_ptr = &*
piolaDiffN.data().begin();
357 t_transformed_diff_n_ptr, &t_transformed_diff_n_ptr[
HVEC0_1],
358 &t_transformed_diff_n_ptr[
HVEC0_2], &t_transformed_diff_n_ptr[
HVEC1_0],
359 &t_transformed_diff_n_ptr[
HVEC1_1], &t_transformed_diff_n_ptr[
HVEC1_2],
360 &t_transformed_diff_n_ptr[
HVEC2_0], &t_transformed_diff_n_ptr[
HVEC2_1],
361 &t_transformed_diff_n_ptr[
HVEC2_2]);
363 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
364 for (
unsigned int bb = 0; bb != nb_base_functions; ++bb) {
370 ++t_transformed_diff_n;
389 const int valid_edges3[] = {1, 1, 1, 0, 0, 0, 0, 0, 0};
390 const int valid_faces3[] = {0, 0, 0, 1, 0, 0, 0, 0, 0};
391 const int valid_edges4[] = {0, 0, 0, 0, 0, 0, 1, 1, 1};
392 const int valid_faces4[] = {0, 0, 0, 0, 1, 0, 0, 0, 0};
394 if (type == MBEDGE) {
395 if (!valid_edges3[side] || valid_edges4[side])
397 }
else if (type == MBTRI) {
398 if (!valid_faces3[side] || valid_faces4[side])
404 for (
unsigned int gg = 0; gg < data.
getN().size1(); ++gg) {
405 for (
int dd = 0; dd < 3; dd++) {
423 if (2 * data.
getN().size2() != data.
getDiffN().size2()) {
427 if (nb_dofs % 3 != 0) {
430 if (nb_dofs > 3 * data.
getN().size2()) {
432 "data inconsistency, side %d type %d", side, type);
434 for (
unsigned int gg = 0; gg < data.
getN().size1(); ++gg) {
435 for (
int dd = 0; dd < 3; dd++) {
436 if ((type == MBTRI && valid_faces3[side]) ||
437 (type == MBEDGE && valid_edges3[side])) {
441 cblas_ddot(nb_dofs / 3, &data.
getDiffN()(gg, 0), 2,
444 cblas_ddot(nb_dofs / 3, &data.
getDiffN()(gg, 1), 2,
446 }
else if ((type == MBTRI && valid_faces4[side]) ||
447 (type == MBEDGE && valid_edges4[side])) {
451 cblas_ddot(nb_dofs / 3, &data.
getDiffN()(gg, 0), 2,
454 cblas_ddot(nb_dofs / 3, &data.
getDiffN()(gg, 1), 2,
476 cblas_dgemv(CblasRowMajor, CblasNoTrans, 3, 3, 1., &*
sPin.data().begin(), 3,
485 cblas_dgemv(CblasRowMajor, CblasNoTrans, 3, 3, 1., &*
sPin.data().begin(), 3,
497 if (moab::CN::Dimension(type) != 2)
502 "Pointer to normal/normals not set");
505 if (normal_is_at_gauss_pts)
508 auto apply_transform_linear_geometry = [&](
auto base,
auto nb_gauss_pts,
509 auto nb_base_functions) {
515 const auto l02 = t_normal(
i) * t_normal(
i);
517 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
518 for (
int bb = 0; bb != nb_base_functions; ++bb) {
519 const auto v = t_base(0);
520 t_base(
i) = (
v / l02) * t_normal(
i);
527 auto apply_transform_nonlinear_geometry = [&](
auto base,
auto nb_gauss_pts,
528 auto nb_base_functions) {
532 &normals_at_pts(0, 0), &normals_at_pts(0, 1), &normals_at_pts(0, 2));
535 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
536 const auto l2 = t_normal(
i) * t_normal(
i);
537 for (
int bb = 0; bb != nb_base_functions; ++bb) {
538 const auto v = t_base(0);
539 t_base(
i) = (
v / l2) * t_normal(
i);
547 if (normal_is_at_gauss_pts) {
551 const auto &base_functions = data.
getN(base);
552 const auto nb_gauss_pts = base_functions.size1();
558 "normalsAtGaussPtsRawPtr has inconsistent number of "
562 const auto nb_base_functions = base_functions.size2() / 3;
563 CHKERR apply_transform_nonlinear_geometry(base, nb_gauss_pts,
571 const auto &base_functions = data.
getN(base);
572 const auto nb_gauss_pts = base_functions.size1();
575 const auto nb_base_functions = base_functions.size2() / 3;
576 CHKERR apply_transform_linear_geometry(base, nb_gauss_pts,
589 const auto type_dim = moab::CN::Dimension(type);
590 if (type_dim != 1 && type_dim != 2)
612 auto &baseN = data.
getN(base);
613 auto &diffBaseN = data.
getDiffN(base);
615 int nb_dofs = baseN.size2() / 3;
616 int nb_gauss_pts = baseN.size1();
619 MatrixDouble diff_piola_n(diffBaseN.size1(), diffBaseN.size2());
621 if (nb_dofs > 0 && nb_gauss_pts > 0) {
632 t_transformed_diff_h_curl(
646 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
649 for (
int ll = 0; ll != nb_dofs; ll++) {
650 t_transformed_h_curl(
i) = t_inv_m(
j,
i) * t_h_curl(
j);
651 t_transformed_diff_h_curl(
i,
k) =
652 t_inv_m(
j,
i) * t_diff_h_curl(
j,
k);
654 ++t_transformed_h_curl;
656 ++t_transformed_diff_h_curl;
662 for (
int gg = 0; gg < nb_gauss_pts; ++gg) {
663 for (
int ll = 0; ll != nb_dofs; ll++) {
664 t_transformed_h_curl(
i) = t_inv_m(
j,
i) * t_h_curl(
j);
665 t_transformed_diff_h_curl(
i,
k) =
666 t_inv_m(
j,
i) * t_diff_h_curl(
j,
k);
668 ++t_transformed_h_curl;
670 ++t_transformed_diff_h_curl;
675 if (cc != nb_gauss_pts * nb_dofs)
679 diffBaseN.swap(diff_piola_n);
695 int nb_gauss_pts = data.
getN().size1();
696 tAngent.resize(nb_gauss_pts, 3,
false);
698 int nb_approx_fun = data.
getN().size2();
699 double *diff = &*data.
getDiffN().data().begin();
703 tAngent.resize(nb_gauss_pts, 3,
false);
707 for (
int dd = 0; dd != 3; dd++) {
708 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
709 tAngent(gg, dd) = cblas_ddot(2, diff, 1, dofs[dd], 3);
716 "Approximated field should be rank 3, i.e. vector in 3d space");
718 for (
int dd = 0; dd != 3; dd++) {
719 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
721 cblas_ddot(nb_dofs / 3, &diff[gg * nb_approx_fun], 1, dofs[dd], 3);
727 "This operator can calculate tangent vector only on edge");
743 const double l0 = t_m(
i) * t_m(
i);
745 auto get_base_at_pts = [&](
auto base) {
752 auto get_tangent_at_pts = [&]() {
759 auto calculate_squared_edge_length = [&]() {
760 std::vector<double> l1;
763 l1.resize(nb_gauss_pts);
764 auto t_m_at_pts = get_tangent_at_pts();
765 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
766 l1[gg] = t_m_at_pts(
i) * t_m_at_pts(
i);
773 auto l1 = calculate_squared_edge_length();
778 const size_t nb_gauss_pts = data.
getN(base).size1();
779 const size_t nb_dofs = data.
getN(base).size2() / 3;
780 if (nb_gauss_pts && nb_dofs) {
781 auto t_h_curl = get_base_at_pts(base);
784 auto t_m_at_pts = get_tangent_at_pts();
785 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
786 const double l0 = l1[gg];
787 for (
int ll = 0; ll != nb_dofs; ll++) {
788 const double val = t_h_curl(0);
789 const double a = val / l0;
790 t_h_curl(
i) = t_m_at_pts(
i) *
a;
797 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
798 for (
int ll = 0; ll != nb_dofs; ll++) {
799 const double val = t_h_curl(0);
800 const double a = val / l0;
801 t_h_curl(
i) = t_m(
i) *
a;
808 if (cc != nb_gauss_pts * nb_dofs)
820 double *ptr = &*data.data().begin();
827OpGetDataAndGradient<3, 3>::getGradAtGaussPtsTensor<3, 3>(
MatrixDouble &data) {
828 double *ptr = &*data.data().begin();
830 &ptr[4], &ptr[5], &ptr[6], &ptr[7],
840 const unsigned int nb_gauss_pts = data.
getN().size1();
841 const unsigned int nb_base_functions = data.
getN().size2();
842 const unsigned int nb_dofs = data.
getFieldData().size();
846 auto t_val = getValAtGaussPtsTensor<3>(dataAtGaussPts);
847 auto t_grad = getGradAtGaussPtsTensor<3, 3>(dataGradAtGaussPts);
850 if (type == MBVERTEX &&
851 data.
getDiffN().data().size() == 3 * nb_base_functions) {
852 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
856 for (; bb != nb_dofs / 3; ++bb) {
857 t_val(
i) += t_data(
i) * t_n;
858 t_grad(
i,
j) += t_data(
i) * t_diff_n(
j);
865 for (; bb != nb_base_functions; ++bb) {
871 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
874 for (; bb != nb_dofs / 3; ++bb) {
875 t_val(
i) += t_data(
i) * t_n;
876 t_grad(
i,
j) += t_data(
i) * t_diff_n(
j);
883 for (; bb != nb_base_functions; ++bb) {
896 const unsigned int nb_gauss_pts = data.
getN().size1();
897 const unsigned int nb_base_functions = data.
getN().size2();
899 const unsigned int nb_dofs = data.
getFieldData().size();
903 double *ptr = &*dataGradAtGaussPts.data().begin();
907 if (type == MBVERTEX &&
908 data.
getDiffN().data().size() == 3 * nb_base_functions) {
909 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
913 for (; bb != nb_dofs / 3; ++bb) {
914 t_val += t_data * t_n;
915 t_grad(
i) += t_data * t_diff_n(
i);
922 for (; bb != nb_base_functions; ++bb) {
928 for (
unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
931 for (; bb != nb_dofs / 3; ++bb) {
932 t_val = t_data * t_n;
933 t_grad(
i) += t_data * t_diff_n(
i);
940 for (; bb != nb_base_functions; ++bb) {
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#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 ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
implementation of Data Operators for Forces and Sources
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
MoFEMErrorCode invertTensor3by3< 3, double, ublas::row_major, DoubleAllocator >(MatrixDouble &jac_data, VectorDouble &det_data, MatrixDouble &inv_jac_data)
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
constexpr IntegrationType I
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
virtual MoFEMErrorCode opLhs(EntitiesFieldData &row_data, EntitiesFieldData &col_data)
virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
virtual MoFEMErrorCode opRhs(EntitiesFieldData &data, const bool error_if_no_base=false)
bool getSymm() const
Get if operator uses symmetry of DOFs or not.
DataOperator(const bool symm=true)
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
virtual std::array< boost::shared_ptr< MatrixDouble >, MaxBernsteinBezierOrder > & getBBDiffNByOrderArray()
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
FieldApproximationBase & getBase()
Get approximation base.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coeffinects.
virtual std::map< std::string, boost::shared_ptr< MatrixDouble > > & getBBDiffNMap()
get hash map of direvarives base function for BB base, key is a field name
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
data structure for finite element entity
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
MatrixDouble & tAngent2_at_GaussPtF3
MatrixDouble & tAngent1_at_GaussPtF3
MatrixDouble & nOrmals_at_GaussPtF3
MatrixDouble & cOords_at_GaussPtF3
MatrixDouble & nOrmals_at_GaussPtF4
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MatrixDouble & cOords_at_GaussPtF4
MoFEMErrorCode calculateNormals()
MatrixDouble & tAngent2_at_GaussPtF4
MatrixDouble & tAngent1_at_GaussPtF4
Get field values and gradients at Gauss points.
MoFEMErrorCode calculateValAndGrad(int side, EntityType type, EntitiesFieldData::EntData &data)
Calculate gradient and values at integration points.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Tensor2< double *, 3, 3 > tInvJac
FTensor::Index< 'j', 3 > j
FTensor::Tensor2< double *, 3, 3 > tInvJac
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Index< 'i', 3 > i
FTensor::Index< 'k', 3 > k
MatrixDouble diffHdivInvJac