v0.9.1
EdgeElementForcesAndSourcesCore.cpp
Go to the documentation of this file.
1 /** \file EdgeElementForcesAndSourcesCoreBase.cpp
2 
3 \brief Implementation of edge element
4 
5 */
6 
7 /* This file is part of MoFEM.
8  * MoFEM is free software: you can redistribute it and/or modify it under
9  * the terms of the GNU Lesser General Public License as published by the
10  * Free Software Foundation, either version 3 of the License, or (at your
11  * option) any later version.
12  *
13  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  * License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
20 
21 #ifdef __cplusplus
22 extern "C" {
23 #endif
24 #include <cblas.h>
25 #include <lapack_wrap.h>
26 // #include <gm_rule.h>
27 #include <quad.h>
28 #ifdef __cplusplus
29 }
30 #endif
31 
32 namespace MoFEM {
33 
35  Interface &m_field)
36  : ForcesAndSourcesCore(m_field),
37  meshPositionsFieldName("MESH_NODE_POSITIONS"),
38  opGetHoTangentOnEdge(tangentAtGaussPts),
39  opCovariantTransform(dIrection, tangentAtGaussPts) {
41  boost::shared_ptr<BaseFunction>(new EdgePolynomialBase());
42 }
43 
47  CHKERR mField.get_moab().get_connectivity(ent, cOnn, numNodes, true);
48  cOords.resize(numNodes * 3, false);
49  CHKERR mField.get_moab().get_coords(cOnn, numNodes, &*cOords.data().begin());
50  dIrection.resize(3, false);
51  cblas_dcopy(3, &cOords[3], 1, &*dIrection.data().begin(), 1);
52  cblas_daxpy(3, -1., &cOords[0], 1, &*dIrection.data().begin(), 1);
53  lEngth = cblas_dnrm2(3, &*dIrection.data().begin(), 1);
55 }
56 
59  int order_data = getMaxDataOrder();
60  int order_row = getMaxRowOrder();
61  int order_col = getMaxColOrder();
62  int rule = getRule(order_row, order_col, order_data);
63 
64  int nb_gauss_pts;
65  if (rule >= 0) {
66  if (rule < QUAD_1D_TABLE_SIZE) {
67  if (QUAD_1D_TABLE[rule]->dim != 1) {
68  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
69  }
70  if (QUAD_1D_TABLE[rule]->order < rule) {
71  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
72  "wrong order %d != %d", QUAD_1D_TABLE[rule]->order, rule);
73  }
74  nb_gauss_pts = QUAD_1D_TABLE[rule]->npoints;
75  gaussPts.resize(2, nb_gauss_pts, false);
76  cblas_dcopy(nb_gauss_pts, &QUAD_1D_TABLE[rule]->points[1], 2,
77  &gaussPts(0, 0), 1);
78  cblas_dcopy(nb_gauss_pts, QUAD_1D_TABLE[rule]->weights, 1,
79  &gaussPts(1, 0), 1);
80  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 2,
81  false);
82  double *shape_ptr =
83  &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
84  cblas_dcopy(2 * nb_gauss_pts, QUAD_1D_TABLE[rule]->points, 1, shape_ptr,
85  1);
86  } else {
87  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
88  "rule > quadrature order %d < %d", rule, QUAD_1D_TABLE_SIZE);
89  nb_gauss_pts = 0;
90  }
91  } else {
92  // If rule is negative, set user defined integration points
93  CHKERR setGaussPts(order_row, order_col, order_data);
94  nb_gauss_pts = gaussPts.size2();
95  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 2,
96  false);
97  if (nb_gauss_pts) {
98  for (int gg = 0; gg != nb_gauss_pts; gg++) {
99  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE)(gg, 0) =
100  N_MBEDGE0(gaussPts(0, gg));
101  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE)(gg, 1) =
102  N_MBEDGE1(gaussPts(0, gg));
103  }
104  }
105  }
107 }
108 
112  const int nb_gauss_pts = gaussPts.size2();
113  coordsAtGaussPts.resize(nb_gauss_pts, 3, false);
114  if (nb_gauss_pts) {
116  &*gaussPts.data().begin());
118  &coordsAtGaussPts(0, 0), &coordsAtGaussPts(0, 1),
119  &coordsAtGaussPts(0, 2));
120  FTensor::Tensor1<double, 3> t_coords_node0(cOords[0], cOords[1], cOords[2]);
121  FTensor::Tensor1<double, 3> t_coords_node1(cOords[3], cOords[4], cOords[5]);
123  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
124  t_coords(i) = N_MBEDGE0(t_ksi) * t_coords_node0(i) +
125  N_MBEDGE1(t_ksi) * t_coords_node1(i);
126  ++t_coords;
127  ++t_ksi;
128  }
129  }
131 }
132 
136 
137  if (dataPtr->get<FieldName_mi_tag>().find(meshPositionsFieldName) !=
138  dataPtr->get<FieldName_mi_tag>().end()) {
142  } else {
143  tangentAtGaussPts.resize(0, 3, false);
144  }
146 }
147 
149  ForcesAndSourcesCore *ptr) {
151  if (!(ptrFE = dynamic_cast<EdgeElementForcesAndSourcesCoreBase *>(ptr)))
152  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
153  "User operator and finite element do not work together");
155 }
156 
157 } // namespace MoFEM
int getMaxRowOrder() const
Get max order of approximation for field in rows.
Deprecated interface functions.
MatrixDouble gaussPts
Matrix of integration points.
static QUAD *const QUAD_1D_TABLE[]
Definition: quad.h:164
virtual moab::Interface & get_moab()=0
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
virtual MoFEMErrorCode opRhs(DataForcesAndSourcesCore &data, const bool error_if_no_base=false)
DataForcesAndSourcesCore & dataH1
void cblas_daxpy(const int N, const double alpha, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_daxpy.c:11
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:506
int npoints
Definition: quad.h:29
void cblas_dcopy(const int N, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_dcopy.c:11
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
#define N_MBEDGE1(x)
edge shape function
Definition: fem_tools.h:80
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513
boost::shared_ptr< const FEDofEntity_multiIndex > dataPtr
Pointer to finite element data dofs.
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Calculate base functions on tetrahedral.
double cblas_dnrm2(const int N, const double *X, const int incX)
Definition: cblas_dnrm2.c:12
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
MoFEMErrorCode getNodesFieldData(DataForcesAndSourcesCore &data, const std::string &field_name) const
Get data on nodes.
#define N_MBEDGE0(x)
edge shape function
Definition: fem_tools.h:79
constexpr int order
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
#define CHKERR
Inline error check.
Definition: definitions.h:601
MoFEMErrorCode getEntityFieldData(DataForcesAndSourcesCore &data, const std::string &field_name, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
MultiIndex Tag for field name.
structure to get information form mofem into DataForcesAndSourcesCore
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
int getMaxColOrder() const
Get max order of approximation for field in columns.
int getMaxDataOrder() const
Get max order of approximation for data fields.
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26
#define QUAD_1D_TABLE_SIZE
Definition: quad.h:163