12 meshPositionsFieldName(
"MESH_NODE_POSITIONS"), aRea(elementMeasure) {}
29 auto get_ftensor_n_diff = [&]() {
37 &
m(0, 0), &
m(0, 1), &
m(0, 2));
42 const size_t nb_gauss_pts =
gaussPts.size2();
58 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
67 }
else if (type == MBQUAD) {
75 const size_t nb_gauss_pts =
gaussPts.size2();
90 auto t_diff = get_ftensor_n_diff();
91 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
92 auto t_coords = get_ftensor_from_vec_3d(
coords);
94 t_t1(
i) += t_coords(
i) * t_diff(N0);
95 t_t2(
i) += t_coords(
i) * t_diff(N1);
107 "Element type not implemented");
124 auto calc_normal = [&](
const double *diff_ptr) {
135 diff_ptr, &diff_ptr[1]);
146 for (
int nn = 0; nn !=
num_nodes; ++nn) {
147 t_t1(
i) += t_coords(
i) * t_diff(N0);
148 t_t2(
i) += t_coords(
i) * t_diff(N1);
153 aRea = sqrt(t_normal(
i) * t_normal(
i));
157 const double *diff_ptr;
161 CHKERR calc_normal(diff_ptr);
167 CHKERR calc_normal(diff_ptr);
171 "Element type not implemented");
186 auto get_rule_by_type = [&]() {
189 return getRule(order_row + 1, order_col + 1, order_data + 1);
191 return getRule(order_row, order_col, order_data);
195 const int rule = get_rule_by_type();
197 auto set_integration_pts_for_tri = [&]() {
208 gaussPts.resize(3, nb_gauss_pts,
false);
209 cblas_dcopy(nb_gauss_pts, &
QUAD_2D_TABLE[rule]->points[1], 3,
211 cblas_dcopy(nb_gauss_pts, &
QUAD_2D_TABLE[rule]->points[2], 3,
219 cblas_dcopy(3 * nb_gauss_pts,
QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
232 auto calc_base_for_tri = [&]() {
234 const size_t nb_gauss_pts =
gaussPts.size2();
237 base.resize(nb_gauss_pts, 3,
false);
238 diff_base.resize(3, 2,
false);
247 auto calc_base_for_quad = [&]() {
249 const size_t nb_gauss_pts =
gaussPts.size2();
252 base.resize(nb_gauss_pts, 4,
false);
253 diff_base.resize(nb_gauss_pts, 8,
false);
254 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
276 CHKERR set_integration_pts_for_tri();
281 CHKERR calc_base_for_quad();
285 "Element type not implemented: %d", type);
291 const size_t nb_gauss_pts =
gaussPts.size2();
295 CHKERR calc_base_for_tri();
298 CHKERR calc_base_for_quad();
302 "Element type not implemented: %d", type);
317 auto dim_type = CN::Dimension(type);
319 auto get_data_on_ents = [&](
auto lower_dim,
auto space) {
322 for (
auto dd = dim_type; dd >= lower_dim; --dd) {
323 int nb_ents = moab::CN::NumSubEntities(type, dd);
324 for (
int ii = 0; ii != nb_ents; ++ii) {
325 auto sub_ent_type = moab::CN::SubEntityType(type, dd, ii);
327 auto &data_on_ent = data->dataOnEntities[sub_ent_type];
330 data->spacesOnEntities[sub_ent_type].set(space);
349 const size_t nb_nodes =
351 double *shape_functions =
353 const size_t nb_gauss_pts =
gaussPts.size2();
355 for (
int gg = 0; gg != nb_gauss_pts; ++gg)
356 for (
int dd = 0; dd != 3; ++dd)
358 nb_nodes, &shape_functions[nb_nodes * gg], 1, &
coords[dd], 3);
368 "User operator and finite element do not work together");
413FaceElementForcesAndSourcesCore::UserDataOperator::loopSideVolumes(
415 return loopSide(fe_name, &fe_method, 3);
423 if (toElePtr->gaussPts.size1() != getGaussPts().size1()) {
425 "Inconsistent numer of weights %ld != %ld",
426 toElePtr->gaussPts.size1(), getGaussPts().size1());
428 if (toElePtr->gaussPts.size2() != getGaussPts().size2()) {
430 "Inconsistent numer of integtaion pts %ld != %ld",
431 toElePtr->gaussPts.size2(), getGaussPts().size2());
436 switch (getFEType()) {
441 "Element type not implemented");
446 &
m(0, 0), &
m(0, 1), &
m(0, 2));
451 auto get_local_coords_triangle = [&]() {
452 std::array<double, 3> ksi0 = {0, 1, 0};
453 std::array<double, 3> ksi1 = {0, 0, 1};
454 std::array<double, 9> ref_shapes;
457 auto &node_coords = getCoords();
458 auto &glob_coords = toElePtr->coords;
459 std::array<double, 6> local_coords;
461 &*node_coords.begin(), &*glob_coords.begin(), 3, local_coords.data());
466 auto get_diff_triangle = [&]() {
469 diff_ptr, &diff_ptr[1]);
473 auto get_jac = [&](
auto &&local_coords,
auto &&t_diff) {
479 for (
int nn = 0; nn != 3; ++nn) {
480 t_jac(
I,
J) += t_local_coords(
I) * t_diff(
J);
488 auto t_mat_tangent = [&](
auto &t1,
auto &t2) {
490 &t1(0), &t1(1), &t1(2), &t2(0), &t2(1), &t2(2)};
494 auto transform = [&](
auto &&t_mat_t,
auto &&t_mat_out_t,
auto &&t_inv_jac) {
500 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
502 t_mat_out_t(
J,
i) = t_mat_t(
I,
i) * t_inv_jac(
I,
J);
509 auto calc_normal = [&](
auto &
n,
auto &t1,
auto &t2) {
513 auto t_t1 = get_ftensor_from_mat_3d(t1);
514 auto t_t2 = get_ftensor_from_mat_3d(t2);
515 auto t_n = get_ftensor_from_mat_3d(
n);
516 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
526 t_mat_tangent(getTangent1AtGaussPts(), getTangent2AtGaussPts()),
527 t_mat_tangent(toElePtr->tangentOneAtGaussPts,
528 toElePtr->tangentTwoAtGaussPts),
529 get_jac(get_local_coords_triangle(), get_diff_triangle())
532 calc_normal(toElePtr->normalsAtGaussPts, toElePtr->tangentOneAtGaussPts,
533 toElePtr->tangentTwoAtGaussPts);
static MoFEMErrorCode get_jac(EntitiesFieldData::EntData &col_data, int gg, MatrixDouble &jac_stress, MatrixDouble &jac)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
#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()
#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< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
FTensor::Index< 'J', DIM1 > J
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
constexpr IntegrationType I
double zeta
Viscous hardening.
#define QUAD_2D_TABLE_SIZE
static QUAD *const QUAD_2D_TABLE[]
FTensor::Index< 'm', 3 > m
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
MatrixDouble tangentOneAtGaussPts
virtual MoFEMErrorCode calculateAreaAndNormal()
Calculate element area and normal of the face.
MoFEMErrorCode operator()()
function is run for every finite element
FaceElementForcesAndSourcesCore(Interface &m_field)
MatrixDouble normalsAtGaussPts
virtual MoFEMErrorCode getSpaceBaseAndOrderOnElement()
Determine approximation space and order of base functions.
virtual MoFEMErrorCode calculateAreaAndNormalAtIntegrationPts()
Calculate element area and normal of the face at integration points.
virtual MoFEMErrorCode setIntegrationPts()
Set integration points.
const EntityHandle * conn
virtual MoFEMErrorCode calculateCoordinatesAtGaussPts()
Calculate coordinate at integration points.
MatrixDouble tangentTwoAtGaussPts
ForcesAndSourcesCore * ptrFE
structure to get information form mofem into EntitiesFieldData
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
int getMaxRowOrder() const
Get max order of approximation for field in rows.
MoFEMErrorCode getSpacesAndBaseOnEntities(EntitiesFieldData &data) const
Get field approximation space and base on entities.
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
EntitiesFieldData & dataH1
MatrixDouble coordsAtGaussPts
coordinated at gauss points
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
MoFEMErrorCode getEntitySense(const EntityType type, boost::ptr_vector< EntitiesFieldData::EntData > &data) const
get sense (orientation) of entity
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > dataOnElement
Entity data on element entity rows fields.
MoFEMErrorCode calHierarchicalBaseFunctionsOnElement()
Calculate base functions.
MatrixDouble gaussPts
Matrix of integration points.
MoFEMErrorCode calBernsteinBezierBaseFunctionsOnElement()
Calculate Bernstein-Bezier base.
int getMaxColOrder() const
Get max order of approximation for field in columns.
EntityType lastEvaluatedElementEntityType
Last evaluated type of element entity.
MoFEMErrorCode getEntityDataOrder(const EntityType type, const FieldSpace space, boost::ptr_vector< EntitiesFieldData::EntData > &data) const
Get the entity data order.
int getMaxDataOrder() const
Get max order of approximation for data fields.
MoFEMErrorCode createDataOnElement(EntityType type)
Create a entity data on element object.
Copy geometry-related data from one element to other.
Calculate base functions on triangle.
Calculate base functions on triangle.
Base volume element used to integrate on skeleton.