v0.14.0
Loading...
Searching...
No Matches
Projection10NodeCoordsOnField.cpp
Go to the documentation of this file.
1/** \file Projection10NodeCoordsOnField.cpp
2 * \brief Project displacements/coordinates from 10 node tetrahedra on
3 * hierarchical approximation base.
4 */
5
6
7
8namespace MoFEM {
9
11 Interface &m_field, std::string field_name, int verb)
12 : mField(m_field), fieldName(field_name), vErbose(verb) {}
13
16 auto field_ptr = mField.get_field_structure(fieldName);
17 if (field_ptr->getApproxBase() == AINSWORTH_BERNSTEIN_BEZIER_BASE) {
18 MOFEM_TAG_AND_LOG("WORLD", Sev::warning, "Projection10NodeCoordsOnField")
19 << "Only working well for first order AINSWORTH_BERNSTEIN_BEZIER_BASE!";
20 MOFEM_LOG_CHANNEL("WORLD");
21 }
23}
24
27
28 if (dofPtr == NULL) {
29 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
30 }
31 if (dofPtr->getName() != fieldName)
33 if (dofPtr->getEntType() == MBVERTEX) {
34 EntityHandle node = dofPtr->getEnt();
35 coords.resize(3);
36 CHKERR mField.get_moab().get_coords(&node, 1, &*coords.data().begin());
37 dofPtr->getFieldData() = coords[dofPtr->getDofCoeffIdx()];
38 if (vErbose > 0) {
39 MOFEM_TAG_AND_LOG_C("SELF", Sev::noisy, "Projection10NodeCoordsOnField",
40 "val = %6.7e\n", dofPtr->getFieldData());
41 MOFEM_LOG_CHANNEL("SELF");
42 }
43
45 }
46 if (dofPtr->getEntType() != MBEDGE) {
48 }
49 if (dofPtr->getEntDofIdx() != dofPtr->getDofCoeffIdx()) {
51 }
52 EntityHandle edge = dofPtr->getEnt();
53 if (type_from_handle(edge) != MBEDGE) {
54 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
55 "this method works only elements which are type of MBEDGE");
56 }
57 // coords
58 int num_nodes;
59 const EntityHandle *conn;
60 CHKERR mField.get_moab().get_connectivity(edge, conn, num_nodes, false);
61 if ((num_nodes != 2) && (num_nodes != 3)) {
62 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
63 "this method works only 4 node and 10 node tets");
64 }
65 coords.resize(num_nodes * 3);
66 CHKERR mField.get_moab().get_coords(conn, num_nodes, &*coords.data().begin());
67 aveMidCoord.resize(3);
68 midNodeCoord.resize(3);
69 for (int dd = 0; dd < 3; dd++) {
70 aveMidCoord[dd] = (coords[0 * 3 + dd] + coords[1 * 3 + dd]) * 0.5;
71 if (num_nodes == 3) {
72 midNodeCoord[dd] = coords[2 * 3 + dd];
73 } else {
74 midNodeCoord[dd] = aveMidCoord[dd];
75 }
76 }
77
78 const auto base = dofPtr->getApproxBase();
79 double edge_shape_function_val;
80 switch (base) {
82 edge_shape_function_val = 0.25;
83 break;
85 edge_shape_function_val = 0.25 * LOBATTO_PHI0(0);
86 break;
88 double L[3];
89 CHKERR Legendre_polynomials(2, 0, NULL, L, NULL, 1);
90 edge_shape_function_val = 0.125 * LOBATTO_PHI0(0);
91 }; break;
92 default:
93 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not yet implemented");
94 }
95
96 diffNodeCoord.resize(3);
97 ublas::noalias(diffNodeCoord) = midNodeCoord - aveMidCoord;
98 dOf.resize(3);
99 ublas::noalias(dOf) = diffNodeCoord / edge_shape_function_val;
100 if (dofPtr->getNbOfCoeffs() > 3) {
101 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
102 "this method works only fields which are rank 3 or lower");
103 }
104 dofPtr->getFieldData() = dOf[dofPtr->getDofCoeffIdx()];
105
107}
108
113
115 std::string _fieldName,
116 bool set_nodes,
117 bool on_coords,
118 std::string on_tag)
119 : Projection10NodeCoordsOnField(m_field, _fieldName), setNodes(set_nodes),
120 onCoords(on_coords), onTag(on_tag), maxApproximationOrder(20) {}
121
122Tag th;
123Field_multiIndex::index<FieldName_mi_tag>::type::iterator field_it;
126
129 if (!onCoords) {
130 if (onTag == "NoNE") {
131 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
132 "tag name not specified");
133 }
135 if (field_it == fieldsPtr->get<FieldName_mi_tag>().end()) {
136 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "field not found %s",
137 fieldName.c_str());
138 }
139 int field_rank = (*field_it)->getNbOfCoeffs();
140 VectorDouble def_VAL = ublas::zero_vector<double>(field_rank);
141 CHKERR mField.get_moab().tag_get_handle(
142 onTag.c_str(), field_rank, MB_TYPE_DOUBLE, th,
143 MB_TAG_CREAT | MB_TAG_SPARSE, &*def_VAL.data().begin());
145 }
146
147 L.resize(maxApproximationOrder + 1);
149 &*L.data().begin(), NULL, 3);
150 K.resize(10);
151 CHKERR LobattoKernel_polynomials(9, 0., NULL, &*K.data().begin(), NULL, 3);
152
154}
155
158 if (dofPtr->getName() != fieldName)
160 if (setNodes) {
161 if (dofPtr->getEntType() == MBVERTEX) {
162 EntityHandle node = dofPtr->getEnt();
163 if (onCoords) {
164 coords.resize(3);
165 CHKERR mField.get_moab().get_coords(&node, 1, &*coords.data().begin());
166 coords[dofPtr->getDofCoeffIdx()] = dofPtr->getFieldData();
167 CHKERR mField.get_moab().set_coords(&node, 1, &*coords.data().begin());
168 } else {
169 int field_rank = (*field_it)->getNbOfCoeffs();
170 if (field_rank != dofPtr->getNbOfCoeffs()) {
171 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
172 "data inconsistency");
173 }
174 double *tag_value;
175 int tag_size;
176 CHKERR mField.get_moab().tag_get_by_ptr(
177 th, &node, 1, (const void **)&tag_value, &tag_size);
178 if (tag_size != dofPtr->getNbOfCoeffs()) {
179 SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
180 }
181 tag_value[dofPtr->getDofCoeffIdx()] = dofPtr->getFieldData();
182 }
183 }
185 }
186 if (dofPtr->getEntType() != MBEDGE) {
188 }
189 EntityHandle edge = dofPtr->getEnt();
190 if (type_from_handle(edge) != MBEDGE) {
191 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
192 "this method works only elements which are type of MBEDGE");
193 }
194
195 int num_nodes;
196 const EntityHandle *conn;
197 CHKERR mField.get_moab().get_connectivity(edge, conn, num_nodes, false);
198 if ((num_nodes != 2) && (num_nodes != 3)) {
199 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
200 "this method works only 4 node and 10 node tets");
201 }
202 if (num_nodes == 2) {
204 }
205
206 if (dofPtr->getDofOrder() >= maxApproximationOrder) {
207 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
208 "too big approximation order, increase constant "
209 "max_ApproximationOrder");
210 }
211 double approx_val = 0;
212 FieldApproximationBase base = dofPtr->getApproxBase();
213 switch (base) {
215 approx_val = 0.25 * L[dofPtr->getDofOrder() - 2] * dofPtr->getFieldData();
216 break;
218 approx_val = 0.25 * K[dofPtr->getDofOrder() - 2] * dofPtr->getFieldData();
219 break;
220 default:
221 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not yet implemented");
222 }
223
224 if (onCoords) {
225 coords.resize(num_nodes * 3);
226 CHKERR mField.get_moab().get_coords(conn, num_nodes,
227 &*coords.data().begin());
228 if (dofPtr->getEntDofIdx() == dofPtr->getDofCoeffIdx()) {
229 // add only one when higher order terms present
230 double ave_mid = (coords[3 * 0 + dofPtr->getDofCoeffIdx()] +
231 coords[3 * 1 + dofPtr->getDofCoeffIdx()]) *
232 0.5;
233 coords[2 * 3 + dofPtr->getDofCoeffIdx()] = ave_mid;
234 }
235 coords[2 * 3 + dofPtr->getDofCoeffIdx()] += approx_val;
236 CHKERR mField.get_moab().set_coords(&conn[2], 1, &coords[3 * 2]);
237 } else {
238 int tag_size;
239 double *tag_value[num_nodes];
240 CHKERR mField.get_moab().tag_get_by_ptr(
241 th, conn, num_nodes, (const void **)tag_value, &tag_size);
242 if (tag_size != dofPtr->getNbOfCoeffs()) {
243 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
244 }
245 if (dofPtr->getEntDofIdx() == dofPtr->getDofCoeffIdx()) {
246 // add only one when higher order terms present
247 double ave_mid = (tag_value[0][dofPtr->getDofCoeffIdx()] +
248 tag_value[1][dofPtr->getDofCoeffIdx()]) *
249 0.5;
250 tag_value[2][dofPtr->getDofCoeffIdx()] = ave_mid;
251 }
252 tag_value[2][dofPtr->getDofCoeffIdx()] += approx_val;
253 }
255}
256}; // namespace MoFEM
#define MOFEM_TAG_AND_LOG_C(channel, severity, tag, format,...)
Tag and log in channel.
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
#define LOBATTO_PHI0(x)
Definitions taken from Hermes2d code.
FieldApproximationBase
approximation base
Definition definitions.h:58
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
@ AINSWORTH_LOBATTO_BASE
Definition definitions.h:62
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition definitions.h:64
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Legendre approximation basis.
PetscErrorCode LobattoKernel_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Kernel Lobatto base functions.
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasVector< double > VectorDouble
Definition Types.hpp:68
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
auto type_from_handle(const EntityHandle h)
get type from entity handle
Field_multiIndex::index< FieldName_mi_tag >::type::iterator field_it
constexpr auto field_name
const Field_multiIndex * fieldsPtr
raw pointer to fields container
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
boost::shared_ptr< DofEntity > dofPtr
MultiIndex Tag for field name.
Projection of edge entities with one mid-node on hierarchical basis.
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode operator()()
function is run for every finite element
Projection10NodeCoordsOnField(Interface &m_field, std::string field_name, int verb=0)
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode preProcess()
function is run at the beginning of loop
ProjectionFieldOn10NodeTet(Interface &m_field, std::string _fieldName, bool set_nodes, bool on_coords, std::string on_tag="NoNE")
MoFEMErrorCode operator()()
function is run for every finite element
Field_multiIndex::index< FieldName_mi_tag >::type::iterator field_it