14 dataH1TrianglesOnly(MBPRISM), dataH1TroughThickness(MBPRISM),
15 opHOCoordsAndNormals(hoCoordsAtGaussPtsF3, nOrmals_at_GaussPtF3,
16 tAngent1_at_GaussPtF3, tAngent2_at_GaussPtF3,
17 hoCoordsAtGaussPtsF4, nOrmals_at_GaussPtF4,
18 tAngent1_at_GaussPtF4, tAngent2_at_GaussPtF4) {
20 "Problem with creation data on element");
29 auto get_fe_coordinates = [&]() {
41 auto calculate_area_of_triangles = [&] {
51 CHKERR get_fe_coordinates();
52 CHKERR calculate_area_of_triangles();
58 auto get_h1_base_data = [&](
auto &
dataH1) {
95 auto set_gauss_points_on_triangle = [&](
int &nb_gauss_pts_on_faces) {
97 int order_triangles_only = 1;
98 int valid_edges[] = {1, 1, 1, 0, 0, 0, 1, 1, 1};
99 for (
unsigned int ee = 0; ee < 9; ee++) {
100 if (!valid_edges[ee])
102 order_triangles_only = std::max(
103 order_triangles_only,
106 for (
unsigned int ff = 3; ff <= 4; ff++) {
107 order_triangles_only = std::max(
108 order_triangles_only,
111 for (
unsigned int qq = 0; qq < 3; qq++) {
112 order_triangles_only = std::max(
113 order_triangles_only,
116 order_triangles_only = std::max(
117 order_triangles_only,
121 nb_gauss_pts_on_faces = 0;
134 cblas_dcopy(nb_gauss_pts_on_faces, &
QUAD_2D_TABLE[rule]->points[1], 3,
136 cblas_dcopy(nb_gauss_pts_on_faces, &
QUAD_2D_TABLE[rule]->points[2], 3,
138 cblas_dcopy(nb_gauss_pts_on_faces,
QUAD_2D_TABLE[rule]->weights, 1,
141 nb_gauss_pts_on_faces, 3,
false);
146 cblas_dcopy(3 * nb_gauss_pts_on_faces,
QUAD_2D_TABLE[rule]->points, 1,
163 if (nb_gauss_pts_on_faces == 0)
166 nb_gauss_pts_on_faces, 3,
false);
167 if (nb_gauss_pts_on_faces) {
188 auto set_gauss_points_through_thickness =
189 [&](
int &nb_gauss_pts_through_thickness) {
191 nb_gauss_pts_through_thickness = 0;
192 int order_thickness = 1;
193 for (
unsigned int ee = 3; ee <= 5; ee++) {
194 order_thickness = std::max(
198 for (
unsigned int qq = 0; qq < 3; qq++) {
199 order_thickness = std::max(
203 order_thickness = std::max(
222 cblas_dcopy(nb_gauss_pts_through_thickness,
225 cblas_dcopy(nb_gauss_pts_through_thickness,
230 "rule > quadrature order %d < %d", rule,
232 nb_gauss_pts_through_thickness = 0;
242 auto set_gauss_points_in_volume = [&](
int nb_gauss_pts_on_faces,
243 int nb_gauss_pts_through_thickness,
246 nb_gauss_pts = nb_gauss_pts_on_faces * nb_gauss_pts_through_thickness;
247 gaussPts.resize(4, nb_gauss_pts,
false);
249 for (
int ggf = 0; ggf < nb_gauss_pts_on_faces; ggf++) {
250 for (
int ggt = 0; ggt < nb_gauss_pts_through_thickness; ggt++, gg++) {
261 int nb_gauss_pts, nb_gauss_pts_through_thickness, nb_gauss_pts_on_faces;
262 CHKERR set_gauss_points_on_triangle(nb_gauss_pts_on_faces);
263 if (!nb_gauss_pts_on_faces)
265 CHKERR set_gauss_points_through_thickness(nb_gauss_pts_through_thickness);
266 if (!nb_gauss_pts_through_thickness)
268 CHKERR set_gauss_points_in_volume(
269 nb_gauss_pts_on_faces, nb_gauss_pts_through_thickness, nb_gauss_pts);
271 auto calc_coordinates_at_triangles = [&]() {
274 for (
int gg = 0; gg < nb_gauss_pts_on_faces; gg++) {
275 for (
int dd = 0; dd < 3; dd++) {
291 auto calc_vertex_base_on_prism = [&]() {
300 for (
int dd = 0; dd != 6; dd++) {
302 for (
int ggf = 0; ggf < nb_gauss_pts_on_faces; ggf++) {
303 int ddd = dd > 2 ? dd - 3 : dd;
312 for (
int ggt = 0; ggt < nb_gauss_pts_through_thickness; ggt++, gg++) {
314 double dzeta, edge_shape;
325 dksi_tri_n * edge_shape;
327 deta_tri_n * edge_shape;
336 auto calc_base_on_prism = [&]() {
347 boost::shared_ptr<BaseFunctionCtx>(
356 "Not yet implemented");
360 "Not yet implemented");
364 "Not yet implemented");
369 "Not yet implemented");
376 auto calc_coordinate_on_prism = [&]() {
380 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
381 for (
int dd = 0; dd < 3; dd++) {
390 auto calculate_volume = [&]() {
395 auto get_t_coords = [&]() {
404 const size_t nb_gauss_pts =
gaussPts.size2();
409 auto t_w = get_t_w();
410 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
412 auto t_coords = get_t_coords();
414 for (
size_t n = 0;
n != 6; ++
n) {
415 t_jac(
i,
j) += t_coords(
i) * t_diff_n(
j);
422 vol += det * t_w / 2;
430 auto calc_ho_triangle_face_normals = [&]() {
433 auto check_field = [&]() {
438 (*field_it)->getId())
454 const auto bit_number =
474 CHKERR calc_coordinates_at_triangles();
475 CHKERR calc_vertex_base_on_prism();
476 CHKERR calc_base_on_prism();
477 CHKERR calc_coordinate_on_prism();
478 vOlume = calculate_volume();
479 CHKERR calc_ho_triangle_face_normals();
492 "User operator and finite element do not work together");
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#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< 'n', SPACE_DIM > n
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Field_multiIndex::index< FieldName_mi_tag >::type::iterator field_it
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
double zeta
Viscous hardening.
#define QUAD_2D_TABLE_SIZE
static QUAD *const QUAD_2D_TABLE[]
static QUAD *const QUAD_1D_TABLE[]
#define QUAD_1D_TABLE_SIZE
const Field_multiIndex * fieldsPtr
raw pointer to fields container
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode opRhs(EntitiesFieldData &data, const bool error_if_no_base=false)
Deprecated interface functions.
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
std::bitset< LASTBASE > bAse
bases on element
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
virtual int getRuleThroughThickness(int order)
MatrixDouble hoCoordsAtGaussPtsF3
MatrixDouble nOrmals_at_GaussPtF3
virtual MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only)
virtual int getRuleTrianglesOnly(int order)
MatrixDouble tAngent2_at_GaussPtF3
EntitiesFieldData dataH1TrianglesOnly
MatrixDouble tAngent2_at_GaussPtF4
MatrixDouble tAngent1_at_GaussPtF3
OpGetCoordsAndNormalsOnPrism opHOCoordsAndNormals
FatPrismElementForcesAndSourcesCore(Interface &m_field)
MatrixDouble gaussPtsTrianglesOnly
virtual MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness)
MatrixDouble hoCoordsAtGaussPtsF4
MatrixDouble gaussPtsThroughThickness
MoFEMErrorCode operator()()
function is run for every finite element
MatrixDouble coordsAtGaussPtsTrianglesOnly
MatrixDouble nOrmals_at_GaussPtF4
MatrixDouble tAngent1_at_GaussPtF4
EntitiesFieldData dataH1TroughThickness
Class used to pass element data to calculate base functions on fat prism.
Calculate base functions on tetrahedralFIXME: Need moab and mofem finite element structure to work (t...
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
MultiIndex Tag for field name.
ForcesAndSourcesCore * ptrFE
structure to get information form mofem into EntitiesFieldData
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
MoFEMErrorCode getSpacesAndBaseOnEntities(EntitiesFieldData &data) const
Get field approximation space and base on entities.
MoFEMErrorCode getEntityFieldData(EntitiesFieldData &data, const int bit_number, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
EntitiesFieldData & dataH1
MatrixDouble coordsAtGaussPts
coordinated at gauss points
MatrixDouble gaussPts
Matrix of integration points.
MoFEMErrorCode getNodesFieldData(EntitiesFieldData &data, const int bit_number) const
Get data on nodes.
MoFEMErrorCode createDataOnElement(EntityType type)
Create a entity data on element object.
MoFEMErrorCode calculateNormals()
Volume finite element base.
std::string meshPositionsFieldName
const EntityHandle * conn