v0.13.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  * MoFEM is free software: you can redistribute it and/or modify it under
7  * the terms of the GNU Lesser General Public License as published by the
8  * Free Software Foundation, either version 3 of the License, or (at your
9  * option) any later version.
10  *
11  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
12  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
14  * License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>
18  */
19 
20 namespace MoFEM {
21 
23  Interface &m_field, std::string field_name, int verb)
24  : mField(m_field), fieldName(field_name), vErbose(verb) {}
25 
28  auto field_ptr = mField.get_field_structure(fieldName);
29  if (field_ptr->getApproxBase() == AINSWORTH_BERNSTEIN_BEZIER_BASE) {
30  MOFEM_TAG_AND_LOG("WORLD", Sev::warning, "Projection10NodeCoordsOnField")
31  << "Only working well for first order AINSWORTH_BERNSTEIN_BEZIER_BASE!";
32  MOFEM_LOG_CHANNEL("WORLD");
33  }
35 }
36 
39 
40  if (dofPtr == NULL) {
41  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
42  }
43  if (dofPtr->getName() != fieldName)
45  if (dofPtr->getEntType() == MBVERTEX) {
46  EntityHandle node = dofPtr->getEnt();
47  coords.resize(3);
48  CHKERR mField.get_moab().get_coords(&node, 1, &*coords.data().begin());
49  dofPtr->getFieldData() = coords[dofPtr->getDofCoeffIdx()];
50  if (vErbose > 0) {
51  MOFEM_TAG_AND_LOG_C("SELF", Sev::noisy, "Projection10NodeCoordsOnField",
52  "val = %6.7e\n", dofPtr->getFieldData());
53  MOFEM_LOG_CHANNEL("SELF");
54  }
55 
57  }
58  if (dofPtr->getEntType() != MBEDGE) {
60  }
61  if (dofPtr->getEntDofIdx() != dofPtr->getDofCoeffIdx()) {
63  }
64  EntityHandle edge = dofPtr->getEnt();
65  if (type_from_handle(edge) != MBEDGE) {
66  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
67  "this method works only elements which are type of MBEDGE");
68  }
69  // coords
70  int num_nodes;
71  const EntityHandle *conn;
72  CHKERR mField.get_moab().get_connectivity(edge, conn, num_nodes, false);
73  if ((num_nodes != 2) && (num_nodes != 3)) {
74  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
75  "this method works only 4 node and 10 node tets");
76  }
77  coords.resize(num_nodes * 3);
78  CHKERR mField.get_moab().get_coords(conn, num_nodes, &*coords.data().begin());
79  aveMidCoord.resize(3);
80  midNodeCoord.resize(3);
81  for (int dd = 0; dd < 3; dd++) {
82  aveMidCoord[dd] = (coords[0 * 3 + dd] + coords[1 * 3 + dd]) * 0.5;
83  if (num_nodes == 3) {
84  midNodeCoord[dd] = coords[2 * 3 + dd];
85  } else {
87  }
88  }
89 
90  const auto base = dofPtr->getApproxBase();
91  double edge_shape_function_val;
92  switch (base) {
94  edge_shape_function_val = 0.25;
95  break;
97  edge_shape_function_val = 0.25 * LOBATTO_PHI0(0);
98  break;
99  case DEMKOWICZ_JACOBI_BASE: {
100  double L[3];
101  CHKERR Legendre_polynomials(2, 0, NULL, L, NULL, 1);
102  edge_shape_function_val = 0.125 * LOBATTO_PHI0(0);
103  }; break;
104  default:
105  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not yet implemented");
106  }
107 
108  diffNodeCoord.resize(3);
109  ublas::noalias(diffNodeCoord) = midNodeCoord - aveMidCoord;
110  dOf.resize(3);
111  ublas::noalias(dOf) = diffNodeCoord / edge_shape_function_val;
112  if (dofPtr->getNbOfCoeffs() != 3) {
113  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
114  "this method works only fields which are rank 3");
115  }
116  dofPtr->getFieldData() = dOf[dofPtr->getDofCoeffIdx()];
117 
119 }
120 
124 }
125 
127  std::string _fieldName,
128  bool set_nodes,
129  bool on_coords,
130  std::string on_tag)
131  : Projection10NodeCoordsOnField(m_field, _fieldName), setNodes(set_nodes),
132  onCoords(on_coords), onTag(on_tag), maxApproximationOrder(20) {}
133 
134 Tag th;
135 Field_multiIndex::index<FieldName_mi_tag>::type::iterator field_it;
138 
141  if (!onCoords) {
142  if (onTag == "NoNE") {
143  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
144  "tag name not specified");
145  }
147  if (field_it == fieldsPtr->get<FieldName_mi_tag>().end()) {
148  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "field not found %s",
149  fieldName.c_str());
150  }
151  int field_rank = (*field_it)->getNbOfCoeffs();
152  VectorDouble def_VAL = ublas::zero_vector<double>(field_rank);
153  CHKERR mField.get_moab().tag_get_handle(
154  onTag.c_str(), field_rank, MB_TYPE_DOUBLE, th,
155  MB_TAG_CREAT | MB_TAG_SPARSE, &*def_VAL.data().begin());
156  MOAB_THROW(rval);
157  }
158 
159  L.resize(maxApproximationOrder + 1);
161  &*L.data().begin(), NULL, 3);
162  K.resize(10);
163  CHKERR LobattoKernel_polynomials(9, 0., NULL, &*K.data().begin(), NULL, 3);
164 
166 }
167 
170  if (dofPtr->getName() != fieldName)
172  if (setNodes) {
173  if (dofPtr->getEntType() == MBVERTEX) {
174  EntityHandle node = dofPtr->getEnt();
175  if (onCoords) {
176  coords.resize(3);
177  CHKERR mField.get_moab().get_coords(&node, 1, &*coords.data().begin());
178  coords[dofPtr->getDofCoeffIdx()] = dofPtr->getFieldData();
179  CHKERR mField.get_moab().set_coords(&node, 1, &*coords.data().begin());
180  } else {
181  int field_rank = (*field_it)->getNbOfCoeffs();
182  if (field_rank != dofPtr->getNbOfCoeffs()) {
183  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
184  "data inconsistency");
185  }
186  double *tag_value;
187  int tag_size;
188  CHKERR mField.get_moab().tag_get_by_ptr(
189  th, &node, 1, (const void **)&tag_value, &tag_size);
190  if (tag_size != dofPtr->getNbOfCoeffs()) {
191  SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
192  }
193  tag_value[dofPtr->getDofCoeffIdx()] = dofPtr->getFieldData();
194  }
195  }
197  }
198  if (dofPtr->getEntType() != MBEDGE) {
200  }
201  EntityHandle edge = dofPtr->getEnt();
202  if (type_from_handle(edge) != MBEDGE) {
203  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
204  "this method works only elements which are type of MBEDGE");
205  }
206 
207  int num_nodes;
208  const EntityHandle *conn;
209  CHKERR mField.get_moab().get_connectivity(edge, conn, num_nodes, false);
210  if ((num_nodes != 2) && (num_nodes != 3)) {
211  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
212  "this method works only 4 node and 10 node tets");
213  }
214  if (num_nodes == 2) {
216  }
217 
218  if (dofPtr->getDofOrder() >= maxApproximationOrder) {
219  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
220  "too big approximation order, increase constant "
221  "max_ApproximationOrder");
222  }
223  double approx_val = 0;
224  FieldApproximationBase base = dofPtr->getApproxBase();
225  switch (base) {
227  approx_val = 0.25 * L[dofPtr->getDofOrder() - 2] * dofPtr->getFieldData();
228  break;
230  approx_val = 0.25 * K[dofPtr->getDofOrder() - 2] * dofPtr->getFieldData();
231  break;
232  default:
233  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not yet implemented");
234  }
235 
236  if (onCoords) {
237  coords.resize(num_nodes * 3);
238  CHKERR mField.get_moab().get_coords(conn, num_nodes,
239  &*coords.data().begin());
240  if (dofPtr->getEntDofIdx() == dofPtr->getDofCoeffIdx()) {
241  // add only one when higher order terms present
242  double ave_mid = (coords[3 * 0 + dofPtr->getDofCoeffIdx()] +
243  coords[3 * 1 + dofPtr->getDofCoeffIdx()]) *
244  0.5;
245  coords[2 * 3 + dofPtr->getDofCoeffIdx()] = ave_mid;
246  }
247  coords[2 * 3 + dofPtr->getDofCoeffIdx()] += approx_val;
248  CHKERR mField.get_moab().set_coords(&conn[2], 1, &coords[3 * 2]);
249  } else {
250  int tag_size;
251  double *tag_value[num_nodes];
252  CHKERR mField.get_moab().tag_get_by_ptr(
253  th, conn, num_nodes, (const void **)tag_value, &tag_size);
254  if (tag_size != dofPtr->getNbOfCoeffs()) {
255  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
256  }
257  if (dofPtr->getEntDofIdx() == dofPtr->getDofCoeffIdx()) {
258  // add only one when higher order terms present
259  double ave_mid = (tag_value[0][dofPtr->getDofCoeffIdx()] +
260  tag_value[1][dofPtr->getDofCoeffIdx()]) *
261  0.5;
262  tag_value[2][dofPtr->getDofCoeffIdx()] = ave_mid;
263  }
264  tag_value[2][dofPtr->getDofCoeffIdx()] += approx_val;
265  }
267 }
268 }; // namespace MoFEM
#define MOFEM_TAG_AND_LOG_C(channel, severity, tag, format,...)
Tag and log in channel.
Definition: LogManager.hpp:363
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:355
#define LOBATTO_PHI0(x)
Definitions taken from Hermes2d code.
FieldApproximationBase
approximation base
Definition: definitions.h:71
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:75
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:79
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:77
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:554
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:45
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
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 Field * get_field_structure(const std::string &name)=0
get field structure
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:287
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
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
auto type_from_handle
get type from entity handle
Definition: Templates.hpp:1441
Field_multiIndex::index< FieldName_mi_tag >::type::iterator field_it
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