v0.9.0
Projection10NodeCoordsOnField.hpp
Go to the documentation of this file.
1 /** \file Projection10NodeCoordsOnField.hpp
2 
3 FIXME: Move code to cpp file.
4 
5 Project displacements/coordinates from 10 node tetrahedra on hierarchical
6 approximation base.
7 
8 This is example how to use MoFEM::DofMethod when some operator for each node
9 need to be applied.
10 
11 */
12 
13 /* This file is part of MoFEM.
14  * MoFEM is free software: you can redistribute it and/or modify it under
15  * the terms of the GNU Lesser General Public License as published by the
16  * Free Software Foundation, either version 3 of the License, or (at your
17  * option) any later version.
18  *
19  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
20  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
22  * License for more details.
23  *
24  * You should have received a copy of the GNU Lesser General Public
25  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
26 
27 #ifndef __PROJECTION10NODECOORDSONFIELD_HPP__
28 #define __PROJECTION10NODECOORDSONFIELD_HPP__
29 
30 using namespace boost::numeric;
31 
32 namespace MoFEM {
33 
34 /** \brief Projection of edge entities with one mid-node on hierarchical basis
35  */
37 
39  std::string fieldName;
40  int vErbose;
41 
42  Projection10NodeCoordsOnField(Interface &m_field, std::string field_name,
43  int verb = 0)
44  : mField(m_field), fieldName(field_name), vErbose(verb) {}
45 
49  }
50 
56 
59  if (dofPtr == NULL) {
60  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
61  }
62  if (dofPtr->getName() != fieldName)
64  if (dofPtr->getEntType() == MBVERTEX) {
65  EntityHandle node = dofPtr->getEnt();
66  coords.resize(3);
67  CHKERR mField.get_moab().get_coords(&node, 1, &*coords.data().begin());
68  dofPtr->getFieldData() = coords[dofPtr->getDofCoeffIdx()];
69  if (vErbose > 0) {
70  PetscPrintf(mField.get_comm(), "val = %6.7e\n", dofPtr->getFieldData());
71  }
73  }
74  if (dofPtr->getEntType() != MBEDGE) {
76  }
77  if (dofPtr->getEntDofIdx() != dofPtr->getDofCoeffIdx()) {
79  }
80  EntityHandle edge = dofPtr->getEnt();
81  if (mField.get_moab().type_from_handle(edge) != MBEDGE) {
82  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
83  "this method works only elements which are type of MBEDGE");
84  }
85  // coords
86  int num_nodes;
87  const EntityHandle *conn;
88  CHKERR mField.get_moab().get_connectivity(edge, conn, num_nodes, false);
89  if ((num_nodes != 2) && (num_nodes != 3)) {
90  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
91  "this method works only 4 node and 10 node tets");
92  }
93  coords.resize(num_nodes * 3);
94  CHKERR mField.get_moab().get_coords(conn, num_nodes,
95  &*coords.data().begin());
96  aveMidCoord.resize(3);
97  midNodeCoord.resize(3);
98  for (int dd = 0; dd < 3; dd++) {
99  aveMidCoord[dd] = (coords[0 * 3 + dd] + coords[1 * 3 + dd]) * 0.5;
100  if (num_nodes == 3) {
101  midNodeCoord[dd] = coords[2 * 3 + dd];
102  } else {
103  midNodeCoord[dd] = aveMidCoord[dd];
104  }
105  }
106  double edge_shape_function_val = 0.25;
107  FieldApproximationBase base = dofPtr->getApproxBase();
108  switch (base) {
110  break;
112  edge_shape_function_val *= LOBATTO_PHI0(0);
113  break;
114  default:
115  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not yet implemented");
116  }
117 
118  diffNodeCoord.resize(3);
119  ublas::noalias(diffNodeCoord) = midNodeCoord - aveMidCoord;
120  dOf.resize(3);
121  ublas::noalias(dOf) = diffNodeCoord / edge_shape_function_val;
122  if (dofPtr->getNbOfCoeffs() != 3) {
123  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
124  "this method works only fields which are rank 3");
125  }
126  dofPtr->getFieldData() = dOf[dofPtr->getDofCoeffIdx()];
128  }
129 
133  }
134 };
135 
137 
138  bool setNodes;
139  bool onCoords;
140  std::string onTag;
141 
143 
144  ProjectionFieldOn10NodeTet(Interface &m_field, std::string _fieldName,
145  bool set_nodes, bool on_coords,
146  std::string on_tag = "NoNE")
147  : Projection10NodeCoordsOnField(m_field, _fieldName), setNodes(set_nodes),
148  onCoords(on_coords), onTag(on_tag), maxApproximationOrder(20) {}
149 
150  Tag th;
151  Field_multiIndex::index<FieldName_mi_tag>::type::iterator field_it;
154 
157  if (!onCoords) {
158  if (onTag == "NoNE") {
159  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
160  "tag name not specified");
161  }
162  field_it = fieldsPtr->get<FieldName_mi_tag>().find(fieldName);
163  if (field_it == fieldsPtr->get<FieldName_mi_tag>().end()) {
164  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
165  "field not found %s", fieldName.c_str());
166  }
167  int field_rank = (*field_it)->getNbOfCoeffs();
168  VectorDouble def_VAL = ublas::zero_vector<double>(field_rank);
169  CHKERR mField.get_moab().tag_get_handle(
170  onTag.c_str(), field_rank, MB_TYPE_DOUBLE, th,
171  MB_TAG_CREAT | MB_TAG_SPARSE, &*def_VAL.data().begin());
172  MOAB_THROW(rval);
173  }
174 
175  L.resize(maxApproximationOrder + 1);
176  CHKERR Legendre_polynomials(maxApproximationOrder, 0., NULL,
177  &*L.data().begin(), NULL, 3);
178  K.resize(10);
179  CHKERR LobattoKernel_polynomials(9, 0., NULL, &*K.data().begin(), NULL, 3);
180 
182  }
183 
186  if (dofPtr->getName() != fieldName)
188  if (setNodes) {
189  if (dofPtr->getEntType() == MBVERTEX) {
190  EntityHandle node = dofPtr->getEnt();
191  if (onCoords) {
192  coords.resize(3);
193  CHKERR mField.get_moab().get_coords(&node, 1,
194  &*coords.data().begin());
195  coords[dofPtr->getDofCoeffIdx()] = dofPtr->getFieldData();
196  CHKERR mField.get_moab().set_coords(&node, 1,
197  &*coords.data().begin());
198  } else {
199  int field_rank = (*field_it)->getNbOfCoeffs();
200  if (field_rank != dofPtr->getNbOfCoeffs()) {
201  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
202  "data inconsistency");
203  }
204  double *tag_value;
205  int tag_size;
206  CHKERR mField.get_moab().tag_get_by_ptr(
207  th, &node, 1, (const void **)&tag_value, &tag_size);
208  if (tag_size != dofPtr->getNbOfCoeffs()) {
209  SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
210  }
211  tag_value[dofPtr->getDofCoeffIdx()] = dofPtr->getFieldData();
212  }
213  }
215  }
216  if (dofPtr->getEntType() != MBEDGE) {
218  }
219  EntityHandle edge = dofPtr->getEnt();
220  if (mField.get_moab().type_from_handle(edge) != MBEDGE) {
221  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
222  "this method works only elements which are type of MBEDGE");
223  }
224 
225  int num_nodes;
226  const EntityHandle *conn;
227  CHKERR mField.get_moab().get_connectivity(edge, conn, num_nodes, false);
228  if ((num_nodes != 2) && (num_nodes != 3)) {
229  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
230  "this method works only 4 node and 10 node tets");
231  }
232  if (num_nodes == 2) {
234  }
235 
236  if (dofPtr->getDofOrder() >= maxApproximationOrder) {
237  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
238  "too big approximation order, increase constant "
239  "max_ApproximationOrder");
240  }
241  double approx_val = 0;
242  FieldApproximationBase base = dofPtr->getApproxBase();
243  switch (base) {
245  approx_val = 0.25 * L[dofPtr->getDofOrder() - 2] * dofPtr->getFieldData();
246  break;
248  approx_val = 0.25 * K[dofPtr->getDofOrder() - 2] * dofPtr->getFieldData();
249  break;
250  default:
251  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not yet implemented");
252  }
253 
254  if (onCoords) {
255  coords.resize(num_nodes * 3);
256  CHKERR mField.get_moab().get_coords(conn, num_nodes,
257  &*coords.data().begin());
258  if (dofPtr->getEntDofIdx() == dofPtr->getDofCoeffIdx()) {
259  // add only one when higher order terms present
260  double ave_mid = (coords[3 * 0 + dofPtr->getDofCoeffIdx()] +
261  coords[3 * 1 + dofPtr->getDofCoeffIdx()]) *
262  0.5;
263  coords[2 * 3 + dofPtr->getDofCoeffIdx()] = ave_mid;
264  }
265  coords[2 * 3 + dofPtr->getDofCoeffIdx()] += approx_val;
266  CHKERR mField.get_moab().set_coords(&conn[2], 1, &coords[3 * 2]);
267  } else {
268  int tag_size;
269  double *tag_value[num_nodes];
270  CHKERR mField.get_moab().tag_get_by_ptr(
271  th, conn, num_nodes, (const void **)tag_value, &tag_size);
272  if (tag_size != dofPtr->getNbOfCoeffs()) {
273  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
274  "data inconsistency");
275  }
276  if (dofPtr->getEntDofIdx() == dofPtr->getDofCoeffIdx()) {
277  // add only one when higher order terms present
278  double ave_mid = (tag_value[0][dofPtr->getDofCoeffIdx()] +
279  tag_value[1][dofPtr->getDofCoeffIdx()]) *
280  0.5;
281  tag_value[2][dofPtr->getDofCoeffIdx()] = ave_mid;
282  }
283  tag_value[2][dofPtr->getDofCoeffIdx()] += approx_val;
284  }
286  }
287 };
288 
289 } // namespace MoFEM
290 
291 #endif // __PROJECTION10NODECOORDSONFIELD_HPP__
PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Legendre approximation basis.
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode preProcess()
function is run at the beginning of loop
virtual moab::Interface & get_moab()=0
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
Data structure to exchange data between mofem and User Loop Methods on entities.It allows to exchange...
#define MOAB_THROW(a)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:602
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Projection of edge entities with one mid-node on hierarchical basis.
Projection10NodeCoordsOnField(Interface &m_field, std::string field_name, int verb=0)
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84
MoFEMErrorCode operator()()
function is run for every finite element
MoFEMErrorCode preProcess()
function is run at the beginning of loop
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
PetscErrorCode LobattoKernel_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Kernel Lobatto base functions.
FieldApproximationBase
approximation base
Definition: definitions.h:144
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:146
ProjectionFieldOn10NodeTet(Interface &m_field, std::string _fieldName, bool set_nodes, bool on_coords, std::string on_tag="NoNE")
Field_multiIndex::index< FieldName_mi_tag >::type::iterator field_it
#define LOBATTO_PHI0(x)
Definitions taken from Hermes2d code.
#define CHKERR
Inline error check.
Definition: definitions.h:596
MoFEMErrorCode operator()()
function is run for every finite element
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:85
MultiIndex Tag for field name.
virtual MPI_Comm & get_comm() const =0
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72