v0.14.0
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 
8 namespace MoFEM {
9 
10 Projection10NodeCoordsOnField::Projection10NodeCoordsOnField(
11  Interface &m_field, std::string field_name, int verb)
12  : mField(m_field), fieldName(field_name), vErbose(verb) {}
13 
14 MoFEMErrorCode Projection10NodeCoordsOnField::preProcess() {
15  MoFEMFunctionBeginHot;
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  }
22  MoFEMFunctionReturnHot(0);
23 }
24 
25 MoFEMErrorCode Projection10NodeCoordsOnField::operator()() {
26  MoFEMFunctionBeginHot;
27 
28  if (dofPtr == NULL) {
29  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
30  }
31  if (dofPtr->getName() != fieldName)
32  MoFEMFunctionReturnHot(0);
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 
44  MoFEMFunctionReturnHot(0);
45  }
46  if (dofPtr->getEntType() != MBEDGE) {
47  MoFEMFunctionReturnHot(0);
48  }
49  if (dofPtr->getEntDofIdx() != dofPtr->getDofCoeffIdx()) {
50  MoFEMFunctionReturnHot(0);
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) {
81  case AINSWORTH_LEGENDRE_BASE:
82  edge_shape_function_val = 0.25;
83  break;
84  case AINSWORTH_LOBATTO_BASE:
85  edge_shape_function_val = 0.25 * LOBATTO_PHI0(0);
86  break;
87  case DEMKOWICZ_JACOBI_BASE: {
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 
106  MoFEMFunctionReturnHot(0);
107 }
108 
109 MoFEMErrorCode Projection10NodeCoordsOnField::postProcess() {
110  MoFEMFunctionBeginHot;
111  MoFEMFunctionReturnHot(0);
112 }
113 
114 ProjectionFieldOn10NodeTet::ProjectionFieldOn10NodeTet(Interface &m_field,
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 
122 Tag th;
123 Field_multiIndex::index<FieldName_mi_tag>::type::iterator field_it;
124 VectorDouble L;
125 VectorDouble K;
126 
127 MoFEMErrorCode ProjectionFieldOn10NodeTet::preProcess() {
128  MoFEMFunctionBeginHot;
129  if (!onCoords) {
130  if (onTag == "NoNE") {
131  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
132  "tag name not specified");
133  }
134  field_it = fieldsPtr->get<FieldName_mi_tag>().find(fieldName);
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());
144  MOAB_THROW(rval);
145  }
146 
147  L.resize(maxApproximationOrder + 1);
148  CHKERR Legendre_polynomials(maxApproximationOrder, 0., NULL,
149  &*L.data().begin(), NULL, 3);
150  K.resize(10);
151  CHKERR LobattoKernel_polynomials(9, 0., NULL, &*K.data().begin(), NULL, 3);
152 
153  MoFEMFunctionReturnHot(0);
154 }
155 
156 MoFEMErrorCode ProjectionFieldOn10NodeTet::operator()() {
157  MoFEMFunctionBeginHot;
158  if (dofPtr->getName() != fieldName)
159  MoFEMFunctionReturnHot(0);
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  }
184  MoFEMFunctionReturnHot(0);
185  }
186  if (dofPtr->getEntType() != MBEDGE) {
187  MoFEMFunctionReturnHot(0);
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) {
203  MoFEMFunctionReturnHot(0);
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) {
214  case AINSWORTH_LEGENDRE_BASE:
215  approx_val = 0.25 * L[dofPtr->getDofOrder() - 2] * dofPtr->getFieldData();
216  break;
217  case AINSWORTH_LOBATTO_BASE:
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  }
254  MoFEMFunctionReturnHot(0);
255 }
256 }; // namespace MoFEM
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::K
VectorDouble K
Definition: Projection10NodeCoordsOnField.cpp:125
MoFEM::Projection10NodeCoordsOnField::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: Projection10NodeCoordsOnField.cpp:109
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
EntityHandle
MoFEM::Projection10NodeCoordsOnField::fieldName
std::string fieldName
Definition: Projection10NodeCoordsOnField.hpp:38
MoFEM::CoreInterface::get_field_structure
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::ProjectionFieldOn10NodeTet::maxApproximationOrder
const int maxApproximationOrder
Definition: Projection10NodeCoordsOnField.hpp:65
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
MoFEM::FieldName_mi_tag
MultiIndex Tag for field name.
Definition: TagMultiIndices.hpp:67
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MOAB_THROW
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:541
MoFEM::Projection10NodeCoordsOnField::mField
Interface & mField
Definition: Projection10NodeCoordsOnField.hpp:37
MoFEM::DofMethod::dofPtr
boost::shared_ptr< DofEntity > dofPtr
Definition: LoopMethods.hpp:505
MoFEM::ProjectionFieldOn10NodeTet::field_it
Field_multiIndex::index< FieldName_mi_tag >::type::iterator field_it
Definition: Projection10NodeCoordsOnField.hpp:71
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::BasicMethod::fieldsPtr
const Field_multiIndex * fieldsPtr
raw pointer to fields container
Definition: LoopMethods.hpp:252
MoFEM::Projection10NodeCoordsOnField::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: Projection10NodeCoordsOnField.cpp:14
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
MoFEM::ProjectionFieldOn10NodeTet::K
VectorDouble K
Definition: Projection10NodeCoordsOnField.hpp:73
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::ProjectionFieldOn10NodeTet::onCoords
bool onCoords
Definition: Projection10NodeCoordsOnField.hpp:62
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::ProjectionFieldOn10NodeTet::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.hpp:72
MoFEM::Projection10NodeCoordsOnField::vErbose
int vErbose
Definition: Projection10NodeCoordsOnField.hpp:39
MoFEM::ProjectionFieldOn10NodeTet::setNodes
bool setNodes
Definition: Projection10NodeCoordsOnField.hpp:61
MoFEM::Projection10NodeCoordsOnField::aveMidCoord
VectorDouble3 aveMidCoord
Definition: Projection10NodeCoordsOnField.hpp:42
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
MoFEM::type_from_handle
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1869
Legendre_polynomials
PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Legendre approximation basis.
Definition: base_functions.c:15
MoFEM::ProjectionFieldOn10NodeTet::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: Projection10NodeCoordsOnField.cpp:127
MoFEM::ProjectionFieldOn10NodeTet::ProjectionFieldOn10NodeTet
ProjectionFieldOn10NodeTet(Interface &m_field, std::string _fieldName, bool set_nodes, bool on_coords, std::string on_tag="NoNE")
Definition: Projection10NodeCoordsOnField.cpp:114
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
LobattoKernel_polynomials
PetscErrorCode LobattoKernel_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Kernel Lobatto base functions.
Definition: base_functions.c:345
LOBATTO_PHI0
#define LOBATTO_PHI0(x)
Definitions taken from Hermes2d code.
Definition: base_functions.h:126
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::Projection10NodeCoordsOnField::dOf
VectorDouble3 dOf
Definition: Projection10NodeCoordsOnField.hpp:45
MoFEM::ProjectionFieldOn10NodeTet::onTag
std::string onTag
Definition: Projection10NodeCoordsOnField.hpp:63
MOFEM_TAG_AND_LOG_C
#define MOFEM_TAG_AND_LOG_C(channel, severity, tag, format,...)
Tag and log in channel.
Definition: LogManager.hpp:370
MoFEM::Projection10NodeCoordsOnField::operator()
MoFEMErrorCode operator()()
function is run for every finite element
Definition: Projection10NodeCoordsOnField.cpp:25
UBlasVector< double >
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
MoFEM::Projection10NodeCoordsOnField::coords
VectorDouble coords
Definition: Projection10NodeCoordsOnField.hpp:41
MoFEM::ProjectionFieldOn10NodeTet::operator()
MoFEMErrorCode operator()()
function is run for every finite element
Definition: Projection10NodeCoordsOnField.cpp:156
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
MoFEM::field_it
Field_multiIndex::index< FieldName_mi_tag >::type::iterator field_it
Definition: Projection10NodeCoordsOnField.cpp:123
MoFEM::ProjectionFieldOn10NodeTet::th
Tag th
Definition: Projection10NodeCoordsOnField.hpp:67
MoFEM::Projection10NodeCoordsOnField::Projection10NodeCoordsOnField
Projection10NodeCoordsOnField(Interface &m_field, std::string field_name, int verb=0)
Definition: Projection10NodeCoordsOnField.cpp:10
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEM::Projection10NodeCoordsOnField::midNodeCoord
VectorDouble3 midNodeCoord
Definition: Projection10NodeCoordsOnField.hpp:43
MoFEM::Projection10NodeCoordsOnField::diffNodeCoord
VectorDouble3 diffNodeCoord
Definition: Projection10NodeCoordsOnField.hpp:44
Generated by Doxygen 1.8.17 and hosted at University of Glasgow