v0.9.1
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
Deprecated interface functions.
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:506
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513
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:607
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:149
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:151
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:601
MoFEMErrorCode operator()()
function is run for every finite element
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:87
MultiIndex Tag for field name.
virtual MPI_Comm & get_comm() const =0
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72