v0.10.0
ContactPrismElementForcesAndSourcesCore.cpp
Go to the documentation of this file.
1 /** \file ContactPrismElementForcesAndSourcesCore.cpp
2 
3 \brief Implementation of contact prism element
4 
5 */
6 
7 /* This file is part of MoFEM.
8  * MoFEM is free software: you can redistribute it and/or modify it under
9  * the terms of the GNU Lesser General Public License as published by the
10  * Free Software Foundation, either version 3 of the License, or (at your
11  * option) any later version.
12  *
13  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  * License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
20 
21 namespace MoFEM {
22 
23 ContactPrismElementForcesAndSourcesCore::
24  ContactPrismElementForcesAndSourcesCore(Interface &m_field)
25  : ForcesAndSourcesCore(m_field),
26  dataOnMaster{
27 
28  nullptr,
29  boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET), // NOFIELD
30  boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET), // H1
31  boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET), // HCURL
32  boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET), // HDIV
33  boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET) // L2
34 
35  },
36  dataOnSlave{
37 
38  nullptr,
39  boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET), // NOFIELD
40  boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET), // H1
41  boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET), // HCURL
42  boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET), // HDIV
43  boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET) // L2
44 
45  },
46  derivedDataOnMaster{
47 
48  nullptr,
49  boost::make_shared<DerivedDataForcesAndSourcesCore>(
50  dataOnMaster[NOFIELD]),
51  boost::make_shared<DerivedDataForcesAndSourcesCore>(dataOnMaster[H1]),
52  boost::make_shared<DerivedDataForcesAndSourcesCore>(
53  dataOnMaster[HCURL]),
54  boost::make_shared<DerivedDataForcesAndSourcesCore>(
55  dataOnMaster[HDIV]),
56  boost::make_shared<DerivedDataForcesAndSourcesCore>(dataOnMaster[L2])
57 
58  },
59  derivedDataOnSlave{
60 
61  nullptr,
62  boost::make_shared<DerivedDataForcesAndSourcesCore>(
63  dataOnSlave[NOFIELD]),
64  boost::make_shared<DerivedDataForcesAndSourcesCore>(dataOnSlave[H1]),
65  boost::make_shared<DerivedDataForcesAndSourcesCore>(
66  dataOnSlave[HCURL]),
67  boost::make_shared<DerivedDataForcesAndSourcesCore>(
68  dataOnSlave[HDIV]),
69  boost::make_shared<DerivedDataForcesAndSourcesCore>(dataOnSlave[L2])
70 
71  },
72  dataH1Master(*dataOnMaster[H1].get()),
73  dataH1Slave(*dataOnSlave[H1].get()),
74  dataNoFieldSlave(*dataOnSlave[NOFIELD].get()),
75  dataNoFieldMaster(*dataOnMaster[NOFIELD].get()),
76  dataHcurlMaster(*dataOnMaster[HCURL].get()),
77  dataHcurlSlave(*dataOnSlave[HCURL].get()),
78  dataHdivMaster(*dataOnMaster[HDIV].get()),
79  dataL2Master(*dataOnMaster[L2].get()),
80  dataHdivSlave(*dataOnSlave[HDIV].get()),
81  dataL2Slave(*dataOnSlave[L2].get()), opContravariantTransform() {
82 
83  getUserPolynomialBase() =
84  boost::shared_ptr<BaseFunction>(new TriPolynomialBase());
85 
86  dataH1Master.dataOnEntities[MBVERTEX].push_back(
87  new DataForcesAndSourcesCore::EntData());
88  dataH1Slave.dataOnEntities[MBVERTEX].push_back(
89  new DataForcesAndSourcesCore::EntData());
90 
91  for (int ee = 0; ee != 3; ++ee) {
92  dataH1Master.dataOnEntities[MBEDGE].push_back(
93  new DataForcesAndSourcesCore::EntData());
94  dataH1Slave.dataOnEntities[MBEDGE].push_back(
95  new DataForcesAndSourcesCore::EntData());
96  }
97 
98  dataH1Master.dataOnEntities[MBTRI].push_back(
99  new DataForcesAndSourcesCore::EntData());
100  dataH1Slave.dataOnEntities[MBTRI].push_back(
101  new DataForcesAndSourcesCore::EntData());
102 
103  dataHdivMaster.dataOnEntities[MBTRI].push_back(
104  new DataForcesAndSourcesCore::EntData());
105  dataHdivSlave.dataOnEntities[MBTRI].push_back(
106  new DataForcesAndSourcesCore::EntData());
107 
108  // Data on elements for proper spaces
109  dataOnMaster[H1]->setElementType(MBTRI);
110  derivedDataOnMaster[H1]->setElementType(MBTRI);
111  dataOnSlave[H1]->setElementType(MBTRI);
112  derivedDataOnSlave[H1]->setElementType(MBTRI);
113 
114  dataOnMaster[HDIV]->setElementType(MBTRI);
115  derivedDataOnMaster[HDIV]->setElementType(MBTRI);
116  dataOnSlave[HDIV]->setElementType(MBTRI);
117  derivedDataOnSlave[HDIV]->setElementType(MBTRI);
118 }
119 
120 MoFEMErrorCode
121 ContactPrismElementForcesAndSourcesCore::setDefaultGaussPts(const int rule) {
122  MoFEMFunctionBegin;
123 
124  if (rule < QUAD_2D_TABLE_SIZE) {
125  if (QUAD_2D_TABLE[rule]->dim != 2) {
126  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
127  }
128  if (QUAD_2D_TABLE[rule]->order < rule) {
129  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
130  "wrong order %d != %d", QUAD_2D_TABLE[rule]->order, rule);
131  }
132  nbGaussPts = QUAD_2D_TABLE[rule]->npoints;
133  // For master and slave
134  gaussPtsMaster.resize(3, nbGaussPts, false);
135  gaussPtsSlave.resize(3, nbGaussPts, false);
136 
137  cblas_dcopy(nbGaussPts, &QUAD_2D_TABLE[rule]->points[1], 3,
138  &gaussPtsMaster(0, 0), 1);
139  cblas_dcopy(nbGaussPts, &QUAD_2D_TABLE[rule]->points[2], 3,
140  &gaussPtsMaster(1, 0), 1);
141  cblas_dcopy(nbGaussPts, QUAD_2D_TABLE[rule]->weights, 1,
142  &gaussPtsMaster(2, 0), 1);
143 
144  gaussPtsSlave = gaussPtsMaster;
145 
146  dataH1Master.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nbGaussPts, 3,
147  false);
148 
149  dataH1Slave.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nbGaussPts, 3,
150  false);
151 
152  double *shape_ptr_master =
153  &*dataH1Master.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
154  cblas_dcopy(3 * nbGaussPts, QUAD_2D_TABLE[rule]->points, 1,
155  shape_ptr_master, 1);
156  double *shape_ptr_slave =
157  &*dataH1Slave.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
158  cblas_dcopy(3 * nbGaussPts, QUAD_2D_TABLE[rule]->points, 1, shape_ptr_slave,
159  1);
160  } else {
161  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
162  "rule > quadrature order %d < %d", rule, QUAD_2D_TABLE_SIZE);
163  nbGaussPts = 0;
164  }
165 
166  MoFEMFunctionReturn(0);
167 }
168 
169 MoFEMErrorCode ContactPrismElementForcesAndSourcesCore::operator()() {
170  MoFEMFunctionBegin;
171 
172  if (numeredEntFiniteElementPtr->getEntType() != MBPRISM)
173  MoFEMFunctionReturnHot(0);
174  CHKERR createDataOnElement();
175 
176  DataForcesAndSourcesCore &data_div = *dataOnElement[HDIV];
177  DataForcesAndSourcesCore &data_curl = *dataOnElement[HCURL];
178  DataForcesAndSourcesCore &data_l2 = *dataOnElement[HCURL];
179 
180  EntityHandle ent = numeredEntFiniteElementPtr->getEnt();
181  auto get_coord_and_normal = [&]() {
182  MoFEMFunctionBegin;
183  int num_nodes;
184  const EntityHandle *conn;
185  CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
186  coords.resize(num_nodes * 3, false);
187  CHKERR mField.get_moab().get_coords(conn, num_nodes,
188  &*coords.data().begin());
189  normal.resize(6, false);
190 
191  for (auto &v : {&tangentMasterOne, &tangentMasterTwo, &tangentSlaveOne,
192  &tangentSlaveTwo}) {
193  v->resize(3);
194  v->clear();
195  }
196 
197  CHKERR Tools::getTriNormal(&coords[0], &normal[0]);
198  CHKERR Tools::getTriNormal(&coords[9], &normal[3]);
199 
200  auto get_vec_ptr = [](VectorDouble &vec_double, int r = 0) {
201  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(
202  &vec_double(r + 0), &vec_double(r + 1), &vec_double(r + 2));
203  };
204 
205  auto t_coords_master = get_vec_ptr(coords);
206  auto t_coords_slave = get_vec_ptr(coords, 9);
207  auto t_normal_master = get_vec_ptr(normal);
208  auto t_normal_slave = get_vec_ptr(normal, 3);
209 
210  auto t_t1_master = get_vec_ptr(tangentMasterOne);
211  auto t_t2_master = get_vec_ptr(tangentMasterTwo);
212  auto t_t1_slave = get_vec_ptr(tangentSlaveOne);
213  auto t_t2_slave = get_vec_ptr(tangentSlaveTwo);
214 
215  const double *diff_ptr = Tools::diffShapeFunMBTRI.data();
216  FTensor::Tensor1<FTensor::PackPtr<const double *, 2>, 2> t_diff(
217  &diff_ptr[0], &diff_ptr[1]);
218 
219  FTensor::Index<'i', 3> i;
220  FTensor::Index<'j', 3> j;
221  FTensor::Index<'k', 3> k;
222 
223  FTensor::Number<0> N0;
224 
225  for (int nn = 0; nn != 3; ++nn) {
226  t_t1_master(i) += t_coords_master(i) * t_diff(N0);
227  t_t1_slave(i) += t_coords_slave(i) * t_diff(N0);
228  ++t_coords_master;
229  ++t_coords_slave;
230  ++t_diff;
231  }
232 
233  aRea[0] = sqrt(t_normal_master(i) * t_normal_master(i)) / 2.;
234  aRea[1] = sqrt(t_normal_slave(i) * t_normal_slave(i)) / 2.;
235 
236  t_t2_master(j) =
237  FTensor::levi_civita(i, j, k) * t_normal_master(k) * t_t1_master(i);
238  t_t2_slave(j) =
239  FTensor::levi_civita(i, j, k) * t_normal_slave(k) * t_t1_slave(i);
240 
241  MoFEMFunctionReturn(0);
242  };
243  CHKERR get_coord_and_normal();
244 
245  CHKERR getSpacesAndBaseOnEntities(dataH1);
246 
247  // H1
248  if ((dataH1.spacesOnEntities[MBVERTEX]).test(H1)) {
249  CHKERR getEntitySense<MBEDGE>(dataH1);
250  CHKERR getEntityDataOrder<MBEDGE>(dataH1, H1);
251  CHKERR getEntitySense<MBTRI>(dataH1);
252  CHKERR getEntityDataOrder<MBTRI>(dataH1, H1);
253  }
254 
255  // Hcurl
256  if ((dataH1.spacesOnEntities[MBEDGE]).test(HCURL)) {
257  CHKERR getEntitySense<MBEDGE>(data_curl);
258  CHKERR getEntityDataOrder<MBEDGE>(data_curl, HCURL);
259  CHKERR getEntitySense<MBTRI>(data_curl);
260  CHKERR getEntityDataOrder<MBTRI>(data_curl, HCURL);
261  }
262 
263  // Hdiv
264  if ((dataH1.spacesOnEntities[MBTRI]).test(HDIV)) {
265  CHKERR getEntitySense<MBTRI>(data_div);
266  CHKERR getEntityDataOrder<MBTRI>(data_div, HDIV);
267  dataHdivSlave.spacesOnEntities[MBTRI].set(HDIV);
268  dataHdivMaster.spacesOnEntities[MBTRI].set(HDIV);
269  data_div.spacesOnEntities[MBTRI].set(HDIV);
270  }
271 
272  // L2
273  if ((dataH1.spacesOnEntities[MBTRI]).test(L2)) {
274  CHKERR getEntitySense<MBTRI>(data_l2);
275  CHKERR getEntityDataOrder<MBTRI>(data_l2, L2);
276  }
277 
278  auto clean_data = [](DataForcesAndSourcesCore &data) {
279  MoFEMFunctionBegin;
280  data.sPace.reset();
281  data.bAse.reset();
282  for (EntityType t = MBVERTEX; t != MBMAXTYPE; ++t) {
283  data.spacesOnEntities[t].reset();
284  data.basesOnEntities[t].reset();
285  }
286  for (int s = 0; s != LASTSPACE; ++s)
287  data.basesOnSpaces[s].reset();
288 
289  MoFEMFunctionReturn(0);
290  };
291 
292  auto copy_data = [](DataForcesAndSourcesCore &data,
293  DataForcesAndSourcesCore &copy_data, const int shift) {
294  MoFEMFunctionBegin;
295 
296  if (shift != 0 && shift != 6) {
297  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
298  "Wrong shift for contact prism element");
299  }
300 
301  data.sPace = copy_data.sPace;
302  data.bAse = copy_data.bAse;
303  for (auto t : {MBVERTEX, MBEDGE, MBTRI}) {
304  data.spacesOnEntities[t] = copy_data.spacesOnEntities[t];
305  data.basesOnEntities[t] = copy_data.basesOnEntities[t];
306  data.basesOnSpaces[t] = copy_data.basesOnSpaces[t];
307  }
308 
309  for (int ii = 0; ii != 3; ++ii) {
310  data.dataOnEntities[MBEDGE][ii].getSense() =
311  copy_data.dataOnEntities[MBEDGE][ii + shift].getSense();
312  data.dataOnEntities[MBEDGE][ii].getDataOrder() =
313  copy_data.dataOnEntities[MBEDGE][ii + shift].getDataOrder();
314  }
315 
316  if (shift == 0) {
317  data.dataOnEntities[MBTRI][0].getSense() =
318  copy_data.dataOnEntities[MBTRI][3].getSense();
319  data.dataOnEntities[MBTRI][0].getDataOrder() =
320  copy_data.dataOnEntities[MBTRI][3].getDataOrder();
321  } else {
322  data.dataOnEntities[MBTRI][0].getSense() =
323  copy_data.dataOnEntities[MBTRI][4].getSense();
324  data.dataOnEntities[MBTRI][0].getDataOrder() =
325  copy_data.dataOnEntities[MBTRI][4].getDataOrder();
326  }
327 
328  MoFEMFunctionReturn(0);
329  };
330 
331  auto copy_data_hdiv = [](DataForcesAndSourcesCore &data,
332  DataForcesAndSourcesCore &copy_data,
333  const int shift) {
334  MoFEMFunctionBegin;
335 
336  if (shift != 3 && shift != 4) {
337  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
338  "Wrong shift for contact prism element");
339  }
340 
341  data.sPace = copy_data.sPace;
342  data.bAse = copy_data.bAse;
343 
344  for (auto t : {MBVERTEX, MBTRI}) {
345  data.spacesOnEntities[t] = copy_data.spacesOnEntities[MBVERTEX];
346  data.basesOnEntities[t] = copy_data.basesOnEntities[MBVERTEX];
347  data.basesOnSpaces[t] = copy_data.basesOnSpaces[MBVERTEX];
348  }
349 
350  auto &cpy_ent_dat = copy_data.dataOnEntities[MBTRI][shift];
351  auto &ent_dat = data.dataOnEntities[MBTRI][0];
352  ent_dat.getBase() = cpy_ent_dat.getBase();
353  ent_dat.getSpace() = cpy_ent_dat.getSpace();
354  ent_dat.getSense() = ent_dat.getSense();
355  ent_dat.getDataOrder() = cpy_ent_dat.getDataOrder();
356 
357  MoFEMFunctionReturn(0);
358  };
359 
360  CHKERR clean_data(dataH1Slave);
361  CHKERR copy_data(dataH1Slave, dataH1, 6);
362  CHKERR clean_data(dataH1Master);
363  CHKERR copy_data(dataH1Master, dataH1, 0);
364 
365  int order_data = getMaxDataOrder();
366  int order_row = getMaxRowOrder();
367  int order_col = getMaxColOrder(); // maybe two different rules?
368  int rule = getRule(order_row, order_col, order_data);
369 
370  if (rule >= 0) {
371 
372  CHKERR setDefaultGaussPts(rule);
373 
374  } else {
375 
376  // Master-Slave
377  if (gaussPtsMaster.size2() != gaussPtsSlave.size2())
378  SETERRQ(
379  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
380  "Number of Gauss Points at Master triangle is different than slave");
381 
382  CHKERR setGaussPts(order_row, order_col, order_data);
383  nbGaussPts = gaussPtsMaster.size2();
384  dataH1Master.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nbGaussPts, 3,
385  false);
386  dataH1Slave.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nbGaussPts, 3,
387  false);
388 
389  if (nbGaussPts) {
390  CHKERR Tools::shapeFunMBTRI<1>(&*dataH1Master.dataOnEntities[MBVERTEX][0]
391  .getN(NOBASE)
392  .data()
393  .begin(),
394  &gaussPtsMaster(0, 0),
395  &gaussPtsMaster(1, 0), nbGaussPts);
396 
397  CHKERR Tools::shapeFunMBTRI<1>(
398  &*dataH1Slave.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin(),
399  &gaussPtsSlave(0, 0), &gaussPtsSlave(1, 0), nbGaussPts);
400  }
401  }
402 
403  if (nbGaussPts == 0)
404  MoFEMFunctionReturnHot(0);
405 
406  // Get coordinates on slave and master
407  {
408  coordsAtGaussPtsMaster.resize(nbGaussPts, 3, false);
409  coordsAtGaussPtsSlave.resize(nbGaussPts, 3, false);
410  for (int gg = 0; gg < nbGaussPts; gg++) {
411  for (int dd = 0; dd < 3; dd++) {
412  coordsAtGaussPtsMaster(gg, dd) = cblas_ddot(
413  3, &dataH1Master.dataOnEntities[MBVERTEX][0].getN(NOBASE)(gg, 0), 1,
414  &coords[dd], 3);
415  coordsAtGaussPtsSlave(gg, dd) = cblas_ddot(
416  3, &dataH1Slave.dataOnEntities[MBVERTEX][0].getN(NOBASE)(gg, 0), 1,
417  &coords[9 + dd], 3);
418  }
419  }
420  }
421 
422  for (int space = HCURL; space != LASTSPACE; ++space)
423  if (dataOnElement[space]) {
424 
425  dataH1Master.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE) =
426  dataOnMaster[H1]->dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE);
427  dataH1Slave.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE) =
428  dataOnSlave[H1]->dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE);
429 
430  if (space == HDIV) {
431 
432  if (dataH1.spacesOnEntities[MBTRI].test(HDIV)) {
433 
434  CHKERR clean_data(dataHdivSlave);
435  CHKERR copy_data_hdiv(dataHdivSlave, data_div, 4);
436  CHKERR clean_data(dataHdivMaster);
437  CHKERR copy_data_hdiv(dataHdivMaster, data_div, 3);
438 
439  dataHdivMaster.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE) =
440  dataOnMaster[H1]->dataOnEntities[MBVERTEX][0].getNSharedPtr(
441  NOBASE);
442  dataHdivSlave.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE) =
443  dataOnSlave[H1]->dataOnEntities[MBVERTEX][0].getNSharedPtr(
444  NOBASE);
445  }
446  }
447  }
448 
449  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
450  if (dataH1.bAse.test(b)) {
451  switch (static_cast<FieldApproximationBase>(b)) {
452  case AINSWORTH_LEGENDRE_BASE:
453  case AINSWORTH_LOBATTO_BASE:
454  if (dataH1.spacesOnEntities[MBVERTEX].test(H1)) {
455  CHKERR getUserPolynomialBase()->getValue(
456  gaussPtsMaster,
457  boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
458  dataH1Master, H1, static_cast<FieldApproximationBase>(b),
459  NOBASE)));
460 
461  CHKERR getUserPolynomialBase()->getValue(
462  gaussPtsSlave,
463  boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
464  dataH1Slave, H1, static_cast<FieldApproximationBase>(b),
465  NOBASE)));
466  }
467  break;
468  case DEMKOWICZ_JACOBI_BASE:
469  if (dataH1.spacesOnEntities[MBTRI].test(HDIV)) {
470 
471  auto get_ftensor_from_mat_3d = [](MatrixDouble &m) {
472  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(
473  &m(0, 0), &m(0, 1), &m(0, 2));
474  };
475 
476  auto get_tensor_vec = [](VectorDouble &n, const int r = 0) {
477  return FTensor::Tensor1<double *, 3>(&n(r + 0), &n(r + 1),
478  &n(r + 2));
479  };
480 
481  CHKERR getUserPolynomialBase()->getValue(
482  gaussPtsMaster,
483  boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
484  dataHdivMaster, HDIV, static_cast<FieldApproximationBase>(b),
485  NOBASE)));
486  CHKERR getUserPolynomialBase()->getValue(
487  gaussPtsSlave,
488  boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
489  dataHdivSlave, HDIV, static_cast<FieldApproximationBase>(b),
490  NOBASE)));
491 
492  opContravariantTransform.normalRawPtr = &normal;
493  opContravariantTransform.normalShift = 0;
494  CHKERR opContravariantTransform.opRhs(dataHdivMaster);
495  opContravariantTransform.normalShift = 3;
496  CHKERR opContravariantTransform.opRhs(dataHdivSlave);
497 
498  }
499 
500  if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
501  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
502  "Not yet implemented");
503  }
504  if (dataH1.spacesOnEntities[MBTET].test(L2)) {
505  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
506  "Not yet implemented");
507  }
508  break;
509  default:
510  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
511  "Not yet implemented");
512  }
513  }
514  }
515 
516  // Iterate over operators
517  CHKERR loopOverOperators();
518 
519  MoFEMFunctionReturn(0);
520 }
521 
522 MoFEMErrorCode ContactPrismElementForcesAndSourcesCore::loopOverOperators() {
523  MoFEMFunctionBegin;
524 
525  const EntityType type = numeredEntFiniteElementPtr->getEntType();
526  constexpr std::array<UserDataOperator::OpType, 2> types{
527  UserDataOperator::OPROW, UserDataOperator::OPCOL};
528  std::array<std::string, 2> last_eval_field_name{std::string(), std::string()};
529 
530  auto oit = opPtrVector.begin();
531  auto hi_oit = opPtrVector.end();
532 
533  for (; oit != hi_oit; oit++) {
534 
535  oit->setPtrFE(this);
536 
537  if (oit->opType == UserDataOperator::OPLAST) {
538 
539  // Set field
540  switch (oit->sPace) {
541  case NOSPACE:
542  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Unknown space");
543  case NOFIELD:
544  case H1:
545  case HCURL:
546  case HDIV:
547  case L2:
548  break;
549  default:
550  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
551  "Not implemented for this space", oit->sPace);
552  }
553 
554  // Reseat all data which all field dependent
555  dataOnMaster[oit->sPace]->resetFieldDependentData();
556  dataOnSlave[oit->sPace]->resetFieldDependentData();
557 
558  // Run operator
559  try {
560  CHKERR oit->opRhs(*dataOnMaster[oit->sPace], false);
561  CHKERR oit->opRhs(*dataOnSlave[oit->sPace], false);
562  }
563  CATCH_OP_ERRORS(*oit);
564 
565  } else {
566  boost::shared_ptr<DataForcesAndSourcesCore> op_master_data[2];
567  boost::shared_ptr<DataForcesAndSourcesCore> op_slave_data[2];
568 
569  for (int ss = 0; ss != 2; ss++) {
570 
571  const std::string field_name =
572  !ss ? oit->rowFieldName : oit->colFieldName;
573  const Field *field_struture = mField.get_field_structure(field_name);
574  const BitFieldId data_id = field_struture->getId();
575  const FieldSpace space = field_struture->getSpace();
576  op_master_data[ss] =
577  !ss ? dataOnMaster[space] : derivedDataOnMaster[space];
578  op_slave_data[ss] =
579  !ss ? dataOnSlave[space] : derivedDataOnSlave[space];
580 
581  if ((oit->getNumeredEntFiniteElementPtr()->getBitFieldIdData() &
582  data_id)
583  .none())
584  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
585  "no data field < %s > on finite element < %s >",
586  field_name.c_str(), feName.c_str());
587 
588  if (oit->getOpType() & types[ss] ||
589  oit->getOpType() & UserDataOperator::OPROWCOL) {
590 
591  switch (space) {
592  case NOSPACE:
593  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown space");
594  break;
595  case NOFIELD:
596  case H1:
597  case HCURL:
598  case HDIV:
599  case L2:
600  break;
601  default:
602  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
603  "Not implemented for this space", space);
604  }
605 
606  if (last_eval_field_name[ss] != field_name) {
607  CHKERR getEntityFieldData(*op_master_data[ss], *op_slave_data[ss],
608  field_name, MBEDGE);
609 
610  if (!ss) {
611  struct Extractor {
612  boost::weak_ptr<EntityCacheNumeredDofs>
613  operator()(boost::shared_ptr<FieldEntity> &e) {
614  return e->entityCacheRowDofs;
615  }
616  };
617 
618  CHKERR getEntityIndices(*op_master_data[ss], *op_slave_data[ss],
619  field_name, getRowFieldEnts(), MBEDGE,
620  MBPRISM, Extractor());
621  } else {
622  struct Extractor {
623  boost::weak_ptr<EntityCacheNumeredDofs>
624  operator()(boost::shared_ptr<FieldEntity> &e) {
625  return e->entityCacheColDofs;
626  }
627  };
628  CHKERR getEntityIndices(*op_master_data[ss], *op_slave_data[ss],
629  field_name, getRowFieldEnts(), MBEDGE,
630  MBPRISM, Extractor());
631  }
632 
633  switch (space) {
634  case NOSPACE:
635  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
636  "unknown space");
637  break;
638  case H1: {
639 
640  auto get_indices = [&](auto &master, auto &slave, auto &ents,
641  auto &&ex) {
642  return getNodesIndices(
643  field_name, ents,
644  master.dataOnEntities[MBVERTEX][0].getIndices(),
645  master.dataOnEntities[MBVERTEX][0].getLocalIndices(),
646  slave.dataOnEntities[MBVERTEX][0].getIndices(),
647  slave.dataOnEntities[MBVERTEX][0].getLocalIndices(), ex);
648  };
649 
650  auto get_data = [&](DataForcesAndSourcesCore &master_data,
651  DataForcesAndSourcesCore &slave_data) {
652  return getNodesFieldData(
653  field_name,
654  master_data.dataOnEntities[MBVERTEX][0].getFieldData(),
655  slave_data.dataOnEntities[MBVERTEX][0].getFieldData(),
656  master_data.dataOnEntities[MBVERTEX][0].getFieldDofs(),
657  slave_data.dataOnEntities[MBVERTEX][0].getFieldDofs(),
658  master_data.dataOnEntities[MBVERTEX][0].getFieldEntities(),
659  slave_data.dataOnEntities[MBVERTEX][0].getFieldEntities(),
660  master_data.dataOnEntities[MBVERTEX][0].getSpace(),
661  slave_data.dataOnEntities[MBVERTEX][0].getSpace(),
662  master_data.dataOnEntities[MBVERTEX][0].getBase(),
663  slave_data.dataOnEntities[MBVERTEX][0].getBase());
664  };
665 
666  if (!ss) {
667 
668  struct Extractor {
669  boost::weak_ptr<EntityCacheNumeredDofs>
670  operator()(boost::shared_ptr<FieldEntity> &e) {
671  return e->entityCacheRowDofs;
672  }
673  };
674 
675  CHKERR get_indices(*op_master_data[ss], *op_slave_data[ss],
676  getRowFieldEnts(), Extractor());
677  } else {
678  struct Extractor {
679  boost::weak_ptr<EntityCacheNumeredDofs>
680  operator()(boost::shared_ptr<FieldEntity> &e) {
681  return e->entityCacheColDofs;
682  }
683  };
684 
685  CHKERR get_indices(*op_master_data[ss], *op_slave_data[ss],
686  getColFieldEnts(), Extractor());
687  }
688 
689  CHKERR get_data(*op_master_data[ss], *op_slave_data[ss]);
690 
691  } break;
692  case HCURL:
693  case HDIV:
694  case L2:
695  break;
696  case NOFIELD:
697  if (!getNinTheLoop()) {
698  // NOFIELD data are the same for each element, can be
699  // retrieved only once
700  if (!ss) {
701  CHKERR getNoFieldRowIndices(*op_master_data[ss], field_name);
702  CHKERR getNoFieldRowIndices(*op_slave_data[ss], field_name);
703  } else {
704  CHKERR getNoFieldColIndices(*op_master_data[ss], field_name);
705  CHKERR getNoFieldColIndices(*op_slave_data[ss], field_name);
706  }
707  CHKERR getNoFieldFieldData(*op_master_data[ss], field_name);
708  CHKERR getNoFieldFieldData(*op_slave_data[ss], field_name);
709  }
710  break;
711  default:
712  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
713  "Not implemented for this space", space);
714  }
715  last_eval_field_name[ss] = field_name;
716  }
717  }
718  }
719 
720  int type;
721 
722  if (UserDataOperator *cast_oit =
723  dynamic_cast<UserDataOperator *>(&*oit)) {
724  type = cast_oit->getFaceType();
725  if (((oit->getOpType() & UserDataOperator::OPROW) ||
726  (oit->getOpType() & UserDataOperator::OPCOL)) &&
727  ((type & UserDataOperator::FACEMASTERMASTER) ||
728  (type & UserDataOperator::FACEMASTERSLAVE) ||
729  (type & UserDataOperator::FACESLAVEMASTER) ||
730  (type & UserDataOperator::FACESLAVESLAVE))) {
731  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
732  "Wrong combination of FaceType and OpType, OPROW or OPCOL "
733  "combined with face-face OpType");
734  }
735 
736  if ((oit->getOpType() & UserDataOperator::OPROWCOL) &&
737  ((type & UserDataOperator::FACEMASTER) ||
738  (type & UserDataOperator::FACESLAVE))) {
739  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
740  "Wrong combination of FaceType and OpType, OPROWCOL "
741  "combined with face-face OpType");
742  }
743 
744  if (!type) {
745  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
746  "Face type is not set");
747  }
748  } else {
749  type = UserDataOperator::FACEMASTER | UserDataOperator::FACESLAVE |
750  UserDataOperator::FACEMASTERMASTER |
751  UserDataOperator::FACEMASTERSLAVE |
752  UserDataOperator::FACESLAVEMASTER |
753  UserDataOperator::FACESLAVESLAVE;
754  }
755 
756  if (oit->getOpType() & UserDataOperator::OPROW &&
757  (type & UserDataOperator::FACEMASTER)) {
758  try {
759  CHKERR oit->opRhs(*op_master_data[0], false);
760  }
761  CATCH_OP_ERRORS(*oit);
762  }
763 
764  if (oit->getOpType() & UserDataOperator::OPROW &&
765  (type & UserDataOperator::FACESLAVE)) {
766  try {
767  CHKERR oit->opRhs(*op_slave_data[0], false);
768  }
769  CATCH_OP_ERRORS(*oit);
770  }
771 
772  if (oit->getOpType() & UserDataOperator::OPCOL &&
773  (type & UserDataOperator::FACEMASTER)) {
774  try {
775  CHKERR oit->opRhs(*op_master_data[1], false);
776  }
777  CATCH_OP_ERRORS(*oit);
778  }
779 
780  if (oit->getOpType() & UserDataOperator::OPCOL &&
781  (type & UserDataOperator::FACESLAVE)) {
782  try {
783  CHKERR oit->opRhs(*op_slave_data[1], false);
784  }
785  CATCH_OP_ERRORS(*oit);
786  }
787 
788  if (oit->getOpType() & UserDataOperator::OPROWCOL &&
789  (type & UserDataOperator::FACEMASTERMASTER)) {
790  try {
791  CHKERR oit->opLhs(*op_master_data[0], *op_master_data[1]);
792  }
793  CATCH_OP_ERRORS(*oit);
794  }
795 
796  if (oit->getOpType() & UserDataOperator::OPROWCOL &&
797  (type & UserDataOperator::FACEMASTERSLAVE)) {
798  try {
799  CHKERR oit->opLhs(*op_master_data[0], *op_slave_data[1]);
800  }
801  CATCH_OP_ERRORS(*oit);
802  }
803 
804  if (oit->getOpType() & UserDataOperator::OPROWCOL &&
805  (type & UserDataOperator::FACESLAVEMASTER)) {
806  try {
807  CHKERR oit->opLhs(*op_slave_data[0], *op_master_data[1]);
808  }
809  CATCH_OP_ERRORS(*oit);
810  }
811 
812  if (oit->getOpType() & UserDataOperator::OPROWCOL &&
813  (type & UserDataOperator::FACESLAVESLAVE)) {
814  try {
815  CHKERR oit->opLhs(*op_slave_data[0], *op_slave_data[1]);
816  }
817  CATCH_OP_ERRORS(*oit);
818  }
819  }
820  }
821  MoFEMFunctionReturn(0);
822 }
823 
824 MoFEMErrorCode ContactPrismElementForcesAndSourcesCore::getEntityFieldData(
825  DataForcesAndSourcesCore &master_data, DataForcesAndSourcesCore &slave_data,
826  const std::string &field_name, const EntityType type_lo,
827  const EntityType type_hi) const {
828  MoFEMFunctionBegin;
829 
830  auto reset_data = [type_lo, type_hi](auto &data) {
831  for (EntityType t = type_lo; t != type_hi; ++t) {
832  for (auto &dat : data.dataOnEntities[t]) {
833  dat.getDataOrder() = 0;
834  dat.getBase() = NOBASE;
835  dat.getSpace() = NOSPACE;
836  dat.getFieldData().resize(0, false);
837  dat.getFieldDofs().resize(0, false);
838  dat.getFieldEntities().resize(0, false);
839  }
840  }
841  };
842  reset_data(master_data);
843  reset_data(slave_data);
844 
845  auto &field_ents = getDataFieldEnts();
846  auto bit_number = mField.get_field_bit_number(field_name);
847  const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
848  bit_number, get_id_for_min_type(type_lo));
849  auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
850  cmp_uid_lo);
851  if (lo != field_ents.end()) {
852  const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
853  bit_number, get_id_for_max_type(type_hi));
854  auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
855  if (lo != hi) {
856  for (auto it = lo; it != hi; ++it)
857  if (auto e = it->lock()) {
858 
859  auto get_data = [&](auto &data, auto type, auto side) {
860  auto &dat = data.dataOnEntities[type][side];
861  auto &ent_field_dofs = dat.getFieldDofs();
862  auto &ent_field_data = dat.getFieldData();
863  dat.getFieldEntities().resize(1, false);
864  dat.getFieldEntities()[0] = e.get();
865  dat.getBase() = e->getApproxBase();
866  dat.getSpace() = e->getSpace();
867  const int ent_order = e->getMaxOrder();
868  dat.getDataOrder() =
869  dat.getDataOrder() > ent_order ? dat.getDataOrder() : ent_order;
870  const auto dof_ent_field_data = e->getEntFieldData();
871  const int nb_dofs_on_ent = e->getNbDofsOnEnt();
872  ent_field_data.resize(nb_dofs_on_ent, false);
873  noalias(ent_field_data) = e->getEntFieldData();
874  ent_field_dofs.resize(nb_dofs_on_ent, false);
875  std::fill(ent_field_dofs.begin(), ent_field_dofs.end(), nullptr);
876  if (auto cache = e->entityCacheDataDofs.lock()) {
877  for (auto dit = cache->loHi[0]; dit != cache->loHi[1]; ++dit) {
878  ent_field_dofs[(*dit)->getEntDofIdx()] =
879  reinterpret_cast<FEDofEntity *>((*dit).get());
880  }
881  }
882  };
883 
884  const EntityType type = e->getEntType();
885  const int side = e->getSideNumberPtr()->side_number;
886 
887  switch (type) {
888  case MBEDGE:
889 
890  if (side < 3)
891  get_data(master_data, type, side);
892  else if (side > 5)
893  get_data(slave_data, type, side - 6);
894 
895  break;
896  case MBTRI:
897 
898  if (side == 3)
899  get_data(master_data, type, 0);
900  if (side == 4)
901  get_data(slave_data, type, 0);
902 
903  break;
904  default:
905  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
906  "Entity type not implemented (FIXME)");
907  };
908 
909  const int brother_side = e->getSideNumberPtr()->brother_side_number;
910  if (brother_side != -1)
911  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
912  "Case with brother side not implemented (FIXME)");
913  }
914  }
915  }
916 
917  MoFEMFunctionReturn(0);
918 }
919 
920 MoFEMErrorCode ContactPrismElementForcesAndSourcesCore::getNodesFieldData(
921  const std::string field_name, VectorDouble &master_nodes_data,
922  VectorDouble &slave_nodes_data, VectorDofs &master_nodes_dofs,
923  VectorDofs &slave_nodes_dofs,
924 
925  VectorFieldEntities &master_field_entities, VectorFieldEntities &slave_field_entities,
926 
927  FieldSpace &master_space,
928  FieldSpace &slave_space, FieldApproximationBase &master_base,
929  FieldApproximationBase &slave_base) const {
930  MoFEMFunctionBegin;
931 
932  auto set_zero = [](auto &nodes_data, auto &nodes_dofs, auto &field_entities) {
933  nodes_data.resize(0, false);
934  nodes_dofs.resize(0, false);
935  field_entities.resize(0, false);
936  };
937  set_zero(master_nodes_data, master_nodes_dofs, master_field_entities);
938  set_zero(slave_nodes_data, slave_nodes_dofs, slave_field_entities);
939 
940  auto field_it = fieldsPtr->get<FieldName_mi_tag>().find(field_name);
941  if (field_it != fieldsPtr->get<FieldName_mi_tag>().end()) {
942 
943  auto bit_number = (*field_it)->getBitNumber();
944  const int nb_dofs_on_vert = (*field_it)->getNbOfCoeffs();
945  master_space = slave_space = (*field_it)->getSpace();
946  master_base = slave_base = (*field_it)->getApproxBase();
947 
948  auto &field_ents = getDataFieldEnts();
949  const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
950  bit_number, get_id_for_min_type<MBVERTEX>());
951  auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
952  cmp_uid_lo);
953  if (lo != field_ents.end()) {
954  const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
955  bit_number, get_id_for_max_type<MBVERTEX>());
956  auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
957  if (lo != hi) {
958 
959  int nb_dofs = 0;
960  for (auto it = lo; it != hi; ++it) {
961  if (auto e = it->lock()) {
962  nb_dofs += e->getEntFieldData().size();
963  }
964  }
965 
966  if (nb_dofs) {
967 
968  auto init_set = [&](auto &nodes_data, auto &nodes_dofs,
969  auto &field_entities) {
970  constexpr int num_nodes = 3;
971  const int max_nb_dofs = nb_dofs_on_vert * num_nodes;
972  nodes_data.resize(max_nb_dofs, false);
973  nodes_dofs.resize(max_nb_dofs, false);
974  field_entities.resize(num_nodes, false);
975  nodes_data.clear();
976  fill(nodes_dofs.begin(), nodes_dofs.end(), nullptr);
977  fill(field_entities.begin(), field_entities.end(), nullptr);
978  };
979 
980  init_set(master_nodes_data, master_nodes_dofs, master_field_entities);
981  init_set(slave_nodes_data, slave_nodes_dofs, slave_field_entities);
982 
983  for (auto it = lo; it != hi; ++it) {
984  if (auto e = it->lock()) {
985 
986  const auto &sn = e->getSideNumberPtr();
987  int side = sn->side_number;
988 
989  auto set_data = [&](auto &nodes_data, auto &nodes_dofs,
990  auto &field_entities, int side, int pos) {
991  field_entities[side] = e.get();
992  if (auto cache = e->entityCacheDataDofs.lock()) {
993  for (auto dit = cache->loHi[0]; dit != cache->loHi[1];
994  ++dit) {
995  const auto dof_idx = (*dit)->getEntDofIdx();
996  nodes_data[pos + dof_idx] = (*dit)->getFieldData();
997  nodes_dofs[pos + dof_idx] =
998  reinterpret_cast<FEDofEntity *>((*dit).get());
999  }
1000  }
1001  };
1002 
1003  if (side < 3)
1004  set_data(master_nodes_data, master_nodes_dofs,
1005  master_field_entities, side, side * nb_dofs_on_vert);
1006  else
1007  set_data(slave_nodes_data, slave_nodes_dofs,
1008  slave_field_entities, (side - 3),
1009  (side - 3) * nb_dofs_on_vert);
1010 
1011  const int brother_side = sn->brother_side_number;
1012  if (brother_side != -1)
1013  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1014  "Not implemented (FIXME please)");
1015  }
1016  }
1017  }
1018  }
1019  }
1020  }
1021 
1022  MoFEMFunctionReturn(0);
1023 }
1024 
1025 template <typename EXTRACTOR>
1026 MoFEMErrorCode ContactPrismElementForcesAndSourcesCore::getEntityIndices(
1027  DataForcesAndSourcesCore &master_data, DataForcesAndSourcesCore &slave_data,
1028  const std::string &field_name, FieldEntity_vector_view &ents_field,
1029  const EntityType type_lo, const EntityType type_hi,
1030  EXTRACTOR &&extractor) const {
1031  MoFEMFunctionBegin;
1032 
1033  auto clear_data = [type_lo, type_hi](auto &data) {
1034  for (EntityType t = type_lo; t != type_hi; ++t) {
1035  for (auto &d : data.dataOnEntities[t]) {
1036  d.getIndices().resize(0, false);
1037  d.getLocalIndices().resize(0, false);
1038  }
1039  }
1040  };
1041 
1042  clear_data(master_data);
1043  clear_data(slave_data);
1044 
1045  auto bit_number = mField.get_field_bit_number(field_name);
1046  const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
1047  bit_number, get_id_for_min_type(type_lo));
1048  auto lo = std::lower_bound(ents_field.begin(), ents_field.end(), lo_uid,
1049  cmp_uid_lo);
1050  if (lo != ents_field.end()) {
1051  const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
1052  bit_number, get_id_for_max_type(type_hi));
1053  auto hi = std::upper_bound(lo, ents_field.end(), hi_uid, cmp_uid_hi);
1054  if (lo != hi) {
1055 
1056  std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
1057 
1058  for (auto it = lo; it != hi; ++it)
1059  if (auto e = it->lock()) {
1060 
1061  const EntityType type = e->getEntType();
1062  const int side = e->getSideNumberPtr()->side_number;
1063  const int brother_side = e->getSideNumberPtr()->brother_side_number;
1064  if (brother_side != -1)
1065  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1066  "Not implemented case");
1067 
1068  auto get_indices = [&](auto &data, const auto type, const auto side) {
1069  if (auto cache = extractor(e).lock()) {
1070  for (auto dit = cache->loHi[0]; dit != cache->loHi[1]; ++dit) {
1071  auto &dof = (**dit);
1072  auto &dat = data.dataOnEntities[type][side];
1073  auto &ent_field_indices = dat.getIndices();
1074  auto &ent_field_local_indices = dat.getLocalIndices();
1075  if (ent_field_indices.empty()) {
1076  const int nb_dofs_on_ent = dof.getNbDofsOnEnt();
1077  ent_field_indices.resize(nb_dofs_on_ent, false);
1078  ent_field_local_indices.resize(nb_dofs_on_ent, false);
1079  std::fill(ent_field_indices.data().begin(),
1080  ent_field_indices.data().end(), -1);
1081  std::fill(ent_field_local_indices.data().begin(),
1082  ent_field_local_indices.data().end(), -1);
1083  }
1084  const int idx = dof.getEntDofIdx();
1085  ent_field_indices[idx] = dof.getPetscGlobalDofIdx();
1086  ent_field_local_indices[idx] = dof.getPetscLocalDofIdx();
1087  }
1088  }
1089  };
1090 
1091  switch (type) {
1092  case MBEDGE:
1093 
1094  if (side < 3)
1095  get_indices(master_data, type, side);
1096  else if (side > 5)
1097  get_indices(slave_data, type, side - 6);
1098 
1099  break;
1100  case MBTRI:
1101 
1102  if (side == 3)
1103  get_indices(master_data, type, 0);
1104  if (side == 4)
1105  get_indices(slave_data, type, 0);
1106 
1107  break;
1108  default:
1109  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1110  "Entity type not implemented");
1111  }
1112  }
1113  }
1114  }
1115 
1116  MoFEMFunctionReturn(0);
1117 }
1118 
1119 template <typename EXTRACTOR>
1120 MoFEMErrorCode ContactPrismElementForcesAndSourcesCore::getNodesIndices(
1121  const std::string field_name, FieldEntity_vector_view &ents_field,
1122  VectorInt &master_nodes_indices, VectorInt &master_local_nodes_indices,
1123  VectorInt &slave_nodes_indices, VectorInt &slave_local_nodes_indices,
1124  EXTRACTOR &&extractor) const {
1125  MoFEMFunctionBegin;
1126 
1127  master_nodes_indices.resize(0, false);
1128  master_local_nodes_indices.resize(0, false);
1129  slave_nodes_indices.resize(0, false);
1130  slave_local_nodes_indices.resize(0, false);
1131 
1132  auto field_it = fieldsPtr->get<FieldName_mi_tag>().find(field_name);
1133  if (field_it != fieldsPtr->get<FieldName_mi_tag>().end()) {
1134 
1135  auto bit_number = (*field_it)->getBitNumber();
1136  const int nb_dofs_on_vert = (*field_it)->getNbOfCoeffs();
1137 
1138  const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
1139  bit_number, get_id_for_min_type<MBVERTEX>());
1140  auto lo = std::lower_bound(ents_field.begin(), ents_field.end(), lo_uid,
1141  cmp_uid_lo);
1142  if (lo != ents_field.end()) {
1143  const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
1144  bit_number, get_id_for_max_type<MBVERTEX>());
1145  auto hi = std::upper_bound(lo, ents_field.end(), hi_uid, cmp_uid_hi);
1146  if (lo != hi) {
1147 
1148  int nb_dofs = 0;
1149  for (auto it = lo; it != hi; ++it) {
1150  if (auto e = it->lock()) {
1151  if (auto cache = extractor(e).lock()) {
1152  nb_dofs += std::distance(cache->loHi[0], cache->loHi[1]);
1153  }
1154  }
1155  }
1156 
1157  if (nb_dofs) {
1158 
1159  constexpr int num_nodes = 3;
1160  const int max_nb_dofs = nb_dofs_on_vert * num_nodes;
1161 
1162  auto set_vec_size = [&](auto &nodes_indices,
1163  auto &local_nodes_indices) {
1164  nodes_indices.resize(max_nb_dofs, false);
1165  local_nodes_indices.resize(max_nb_dofs, false);
1166  std::fill(nodes_indices.begin(), nodes_indices.end(), -1);
1167  std::fill(local_nodes_indices.begin(), local_nodes_indices.end(),
1168  -1);
1169  };
1170 
1171  set_vec_size(master_nodes_indices, master_local_nodes_indices);
1172  set_vec_size(slave_nodes_indices, slave_local_nodes_indices);
1173 
1174  for (auto it = lo; it != hi; ++it) {
1175  if (auto e = it->lock()) {
1176 
1177  const int side = e->getSideNumberPtr()->side_number;
1178  const int brother_side =
1179  e->getSideNumberPtr()->brother_side_number;
1180  if (brother_side != -1)
1181  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1182  "Not implemented case");
1183 
1184  auto get_indices = [&](auto &nodes_indices,
1185  auto &local_nodes_indices,
1186  const auto side) {
1187  if (auto cache = extractor(e).lock()) {
1188  for (auto dit = cache->loHi[0]; dit != cache->loHi[1];
1189  ++dit) {
1190  const int idx = (*dit)->getPetscGlobalDofIdx();
1191  const int local_idx = (*dit)->getPetscLocalDofIdx();
1192  const int pos =
1193  side * nb_dofs_on_vert + (*dit)->getDofCoeffIdx();
1194  nodes_indices[pos] = idx;
1195  local_nodes_indices[pos] = local_idx;
1196  }
1197  }
1198  };
1199 
1200  if (side < 3)
1201  get_indices(master_nodes_indices, master_local_nodes_indices,
1202  side);
1203  else if (side > 2)
1204  get_indices(slave_nodes_indices, slave_local_nodes_indices,
1205  side - 3);
1206  else
1207  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1208  "Impossible case");
1209  }
1210  }
1211  }
1212  }
1213  }
1214  }
1215 
1216  MoFEMFunctionReturn(0);
1217 }
1218 
1219 } // namespace MoFEM
NOSPACE
@ NOSPACE
Definition: definitions.h:175
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
MoFEM::ContactPrismElementForcesAndSourcesCore::UserDataOperator::FACESLAVEMASTER
@ FACESLAVEMASTER
Definition: ContactPrismElementForcesAndSourcesCore.hpp:93
MoFEM::BasicMethod::getNinTheLoop
int getNinTheLoop() const
get number of evaluated element in the loop
Definition: LoopMethods.hpp:251
MoFEM::DataForcesAndSourcesCore::basesOnEntities
std::array< std::bitset< LASTBASE >, MBMAXTYPE > basesOnEntities
bases on entity types
Definition: DataStructures.hpp:1036
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:26
EntityHandle
MoFEM::ContactPrismElementForcesAndSourcesCore::derivedDataOnMaster
const std::array< boost::shared_ptr< DataForcesAndSourcesCore >, LASTSPACE > derivedDataOnMaster
Entity data on element entity columns fields.
Definition: ContactPrismElementForcesAndSourcesCore.hpp:293
H1
@ H1
continuous field
Definition: definitions.h:177
MoFEM::FEDofEntity
keeps information about indexed dofs for the finite element
Definition: DofsMultiIndices.hpp:281
LASTBASE
@ LASTBASE
Definition: definitions.h:161
MoFEM::ForcesAndSourcesCore::getNoFieldRowIndices
MoFEMErrorCode getNoFieldRowIndices(DataForcesAndSourcesCore &data, const std::string &field_name) const
get col NoField indices
Definition: ForcesAndSourcesCore.cpp:468
MoFEM::ForcesAndSourcesCore::mField
Interface & mField
Definition: ForcesAndSourcesCore.hpp:36
MoFEM::ForcesAndSourcesCore
structure to get information form mofem into DataForcesAndSourcesCore
Definition: ForcesAndSourcesCore.hpp:34
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
MoFEM::VectorFieldEntities
ublas::vector< FieldEntity *, FieldEntAllocator > VectorFieldEntities
Definition: DataStructures.hpp:43
MoFEM::cmp_uid_hi
static auto cmp_uid_hi(const UId &b, const boost::weak_ptr< FieldEntity > &a)
Definition: ForcesAndSourcesCore.cpp:43
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPCOL
@ OPCOL
Definition: ForcesAndSourcesCore.hpp:100
n
static Index< 'n', 3 > n
Definition: BasicFeTools.hpp:78
NOBASE
@ NOBASE
Definition: definitions.h:151
FTensor::d
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
MoFEM::Types::MatrixDouble
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
L2
@ L2
field with C-1 continuity
Definition: definitions.h:180
MoFEM::ForcesAndSourcesCore::setGaussPts
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
Definition: ForcesAndSourcesCore.cpp:1597
MoFEM::ContactPrismElementForcesAndSourcesCore::coords
VectorDouble coords
Definition: ContactPrismElementForcesAndSourcesCore.hpp:259
QUAD_2D_TABLE
static QUAD *const QUAD_2D_TABLE[]
Definition: quad.h:175
MoFEM::ContactPrismElementForcesAndSourcesCore::aRea
std::array< double, 2 > aRea
Array storing master and slave faces areas.
Definition: ContactPrismElementForcesAndSourcesCore.hpp:255
MoFEM::ContactPrismElementForcesAndSourcesCore::ContactPrismElementForcesAndSourcesCore
ContactPrismElementForcesAndSourcesCore(Interface &m_field)
Definition: ContactPrismElementForcesAndSourcesCore.cpp:24
MoFEM::Tools::diffShapeFunMBTRI
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:87
MoFEM::BasicMethod::fieldsPtr
const Field_multiIndex * fieldsPtr
raw pointer to fields container
Definition: LoopMethods.hpp:269
MoFEM::ContactPrismElementForcesAndSourcesCore::dataHdivMaster
DataForcesAndSourcesCore & dataHdivMaster
Definition: ContactPrismElementForcesAndSourcesCore.hpp:304
FTensor::Tensor1
Definition: Tensor1_value.hpp:9
dim
const int dim
Definition: elec_phys_2D.cpp:12
QUAD_::npoints
int npoints
Definition: quad.h:29
MoFEM::ContactPrismElementForcesAndSourcesCore::dataOnMaster
const std::array< boost::shared_ptr< DataForcesAndSourcesCore >, LASTSPACE > dataOnMaster
Entity data on element entity rows fields.
Definition: ContactPrismElementForcesAndSourcesCore.hpp:282
QUAD_2D_TABLE_SIZE
#define QUAD_2D_TABLE_SIZE
Definition: quad.h:174
MoFEM::ContactPrismElementForcesAndSourcesCore::UserDataOperator::FACEMASTER
@ FACEMASTER
Definition: ContactPrismElementForcesAndSourcesCore.hpp:89
MoFEM::ContactPrismElementForcesAndSourcesCore::opContravariantTransform
OpSetContravariantPiolaTransformOnFace opContravariantTransform
Definition: ContactPrismElementForcesAndSourcesCore.hpp:272
MoFEM::cmp_uid_lo
static auto cmp_uid_lo(const boost::weak_ptr< FieldEntity > &a, const UId &b)
Definition: ForcesAndSourcesCore.cpp:32
MoFEM::ContactPrismElementForcesAndSourcesCore::UserDataOperator::FACEMASTERMASTER
@ FACEMASTERMASTER
Definition: ContactPrismElementForcesAndSourcesCore.hpp:91
MoFEM::ContactPrismElementForcesAndSourcesCore::UserDataOperator::FACEMASTERSLAVE
@ FACEMASTERSLAVE
Definition: ContactPrismElementForcesAndSourcesCore.hpp:92
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
FTensor::Number< 0 >
MoFEM::ForcesAndSourcesCore::getMaxRowOrder
int getMaxRowOrder() const
Get max order of approximation for field in rows.
Definition: ForcesAndSourcesCore.cpp:138
N0
static Number< 0 > N0
Definition: BasicFeTools.hpp:89
MoFEM::FEMethod::getColFieldEnts
auto & getColFieldEnts() const
Definition: LoopMethods.hpp:419
EntData
DataForcesAndSourcesCore::EntData EntData
Definition: quad_polynomial_approximation.cpp:27
MoFEM::ForcesAndSourcesCore::getUserPolynomialBase
auto & getUserPolynomialBase()
Get the User Polynomial Base object.
Definition: ForcesAndSourcesCore.hpp:439
MoFEM::ContactPrismElementForcesAndSourcesCore::tangentSlaveTwo
VectorDouble tangentSlaveTwo
Definition: ContactPrismElementForcesAndSourcesCore.hpp:270
MoFEM::DataForcesAndSourcesCore::sPace
std::bitset< LASTSPACE > sPace
spaces on element
Definition: DataStructures.hpp:1029
FieldSpace
FieldSpace
approximation spaces
Definition: definitions.h:174
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPLAST
@ OPLAST
Definition: ForcesAndSourcesCore.hpp:102
m
static Index< 'm', 3 > m
Definition: BasicFeTools.hpp:77
MoFEM::CoreInterface::get_field_structure
virtual const Field * get_field_structure(const std::string &name)=0
get field structure
MoFEM::OpSetContravariantPiolaTransformOnFace::normalShift
int normalShift
Shift in vector for linear geometry.
Definition: DataOperators.hpp:521
MoFEM::ContactPrismElementForcesAndSourcesCore::loopOverOperators
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
Definition: ContactPrismElementForcesAndSourcesCore.cpp:522
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:604
MoFEM::ContactPrismElementForcesAndSourcesCore::normal
VectorDouble normal
vector storing vector normal to master or slave element
Definition: ContactPrismElementForcesAndSourcesCore.hpp:258
MoFEM::ContactPrismElementForcesAndSourcesCore::gaussPtsSlave
MatrixDouble gaussPtsSlave
Definition: ContactPrismElementForcesAndSourcesCore.hpp:267
MoFEM::ContactPrismElementForcesAndSourcesCore::nbGaussPts
int nbGaussPts
Definition: ContactPrismElementForcesAndSourcesCore.hpp:406
MoFEM::ContactPrismElementForcesAndSourcesCore::derivedDataOnSlave
const std::array< boost::shared_ptr< DataForcesAndSourcesCore >, LASTSPACE > derivedDataOnSlave
Definition: ContactPrismElementForcesAndSourcesCore.hpp:295
MoFEM::ForcesAndSourcesCore::dataH1
DataForcesAndSourcesCore & dataH1
Definition: ForcesAndSourcesCore.hpp:807
MoFEM::ForcesAndSourcesCore::getRule
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
Definition: ForcesAndSourcesCore.cpp:1591
convert.type
type
Definition: convert.py:66
MoFEM::ContactPrismElementForcesAndSourcesCore::getNodesFieldData
MoFEMErrorCode getNodesFieldData(const std::string field_name, VectorDouble &master_nodes_data, VectorDouble &slave_nodes_data, VectorDofs &master_nodes_dofs, VectorDofs &slave_nodes_dofs, VectorFieldEntities &master_field_entities, VectorFieldEntities &slave_field_entities, FieldSpace &master_space, FieldSpace &slave_space, FieldApproximationBase &master_base, FieldApproximationBase &slave_base) const
function that gets nodes field data.
Definition: ContactPrismElementForcesAndSourcesCore.cpp:920
MoFEM::FieldEntity::getLocalUniqueIdCalculate
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
Definition: FieldEntsMultiIndices.hpp:144
MoFEM::ForcesAndSourcesCore::getNoFieldColIndices
MoFEMErrorCode getNoFieldColIndices(DataForcesAndSourcesCore &data, const std::string &field_name) const
get col NoField indices
Definition: ForcesAndSourcesCore.cpp:479
FTensor::Index< 'i', 3 >
MoFEM::ContactPrismElementForcesAndSourcesCore::UserDataOperator::FACESLAVE
@ FACESLAVE
Definition: ContactPrismElementForcesAndSourcesCore.hpp:90
MoFEM::FEMethod::getDataFieldEnts
const FieldEntity_vector_view & getDataFieldEnts() const
Definition: LoopMethods.hpp:400
ElectroPhysiology::b
const double b
Definition: ElecPhysOperators.hpp:30
MoFEM::ContactPrismElementForcesAndSourcesCore::getEntityFieldData
MoFEMErrorCode getEntityFieldData(DataForcesAndSourcesCore &master_data, DataForcesAndSourcesCore &slave_data, const std::string &field_name, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
function that gets entity field data.
Definition: ContactPrismElementForcesAndSourcesCore.cpp:824
StdRDOperators::order
const int order
Definition: rd_stdOperators.hpp:28
MoFEM::get_id_for_max_type
EntityHandle get_id_for_max_type()
Definition: RefEntsMultiIndices.hpp:26
LASTSPACE
@ LASTSPACE
FieldSpace in [ 0, LASTSPACE )
Definition: definitions.h:181
CATCH_OP_ERRORS
#define CATCH_OP_ERRORS(OP)
Definition: ForcesAndSourcesCore.cpp:1280
MoFEM::ContactPrismElementForcesAndSourcesCore::operator()
MoFEMErrorCode operator()()
function is run for every finite element
Definition: ContactPrismElementForcesAndSourcesCore.cpp:169
FTensor::levi_civita
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
Definition: Levi_Civita.hpp:617
MoFEM::Field::getSpace
FieldSpace getSpace() const
Get field approximation space.
Definition: FieldMultiIndices.hpp:268
MoFEM::FieldName_mi_tag
MultiIndex Tag for field name.
Definition: TagMultiIndices.hpp:80
MoFEM::ContactPrismElementForcesAndSourcesCore::gaussPtsMaster
MatrixDouble gaussPtsMaster
Definition: ContactPrismElementForcesAndSourcesCore.hpp:265
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:154
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MoFEM::ContactPrismElementForcesAndSourcesCore::coordsAtGaussPtsSlave
MatrixDouble coordsAtGaussPtsSlave
Definition: ContactPrismElementForcesAndSourcesCore.hpp:262
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROWCOL
@ OPROWCOL
Definition: ForcesAndSourcesCore.hpp:101
MoFEM::FEMethod::numeredEntFiniteElementPtr
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
Definition: LoopMethods.hpp:389
MoFEM::ContactPrismElementForcesAndSourcesCore::coordsAtGaussPtsMaster
MatrixDouble coordsAtGaussPtsMaster
Definition: ContactPrismElementForcesAndSourcesCore.hpp:260
MoFEM::ContactPrismElementForcesAndSourcesCore::dataH1Master
DataForcesAndSourcesCore & dataH1Master
Definition: ContactPrismElementForcesAndSourcesCore.hpp:297
MoFEM::ContactPrismElementForcesAndSourcesCore::tangentMasterOne
VectorDouble tangentMasterOne
Definition: ContactPrismElementForcesAndSourcesCore.hpp:271
MoFEM::DataForcesAndSourcesCore
data structure for finite element entity
Definition: DataStructures.hpp:52
MoFEM::DataOperator::opRhs
virtual MoFEMErrorCode opRhs(DataForcesAndSourcesCore &data, const bool error_if_no_base=false)
Definition: DataOperators.cpp:171
MoFEM::Field::getId
const BitFieldId & getId() const
Get unique field id.
Definition: FieldMultiIndices.hpp:246
MoFEM::ForcesAndSourcesCore::getNoFieldFieldData
MoFEMErrorCode getNoFieldFieldData(const std::string field_name, VectorDouble &ent_field_data, VectorDofs &ent_field_dofs, VectorFieldEntities &ent_field) const
Get field data on nodes.
Definition: ForcesAndSourcesCore.cpp:811
MoFEM::CoreInterface::get_field_bit_number
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
MoFEM::OpSetContravariantPiolaTransformOnFace::normalRawPtr
const VectorDouble * normalRawPtr
Definition: DataOperators.hpp:509
MoFEM::ForcesAndSourcesCore::dataOnElement
const std::array< boost::shared_ptr< DataForcesAndSourcesCore >, LASTSPACE > dataOnElement
Entity data on element entity rows fields.
Definition: ForcesAndSourcesCore.hpp:797
MoFEM::Field
Provide data structure for (tensor) field approximation.
Definition: FieldMultiIndices.hpp:64
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
Definition: ForcesAndSourcesCore.hpp:99
MoFEM::ContactPrismElementForcesAndSourcesCore::tangentSlaveOne
VectorDouble tangentSlaveOne
Definition: ContactPrismElementForcesAndSourcesCore.hpp:270
MoFEM::DataForcesAndSourcesCore::basesOnSpaces
std::array< std::bitset< LASTBASE >, LASTSPACE > basesOnSpaces
base on spaces
Definition: DataStructures.hpp:1038
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:158
MoFEM::ContactPrismElementForcesAndSourcesCore::dataH1Slave
DataForcesAndSourcesCore & dataH1Slave
Definition: ContactPrismElementForcesAndSourcesCore.hpp:298
MoFEM::DataForcesAndSourcesCore::dataOnEntities
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Definition: DataStructures.hpp:1040
MoFEM::ContactPrismElementForcesAndSourcesCore::dataOnSlave
const std::array< boost::shared_ptr< DataForcesAndSourcesCore >, LASTSPACE > dataOnSlave
Definition: ContactPrismElementForcesAndSourcesCore.hpp:284
MoFEM::ForcesAndSourcesCore::getMaxDataOrder
int getMaxDataOrder() const
Get max order of approximation for data fields.
Definition: ForcesAndSourcesCore.cpp:127
MoFEM::get_id_for_min_type
EntityHandle get_id_for_min_type()
Definition: RefEntsMultiIndices.hpp:31
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::ForcesAndSourcesCore::getSpacesAndBaseOnEntities
MoFEMErrorCode getSpacesAndBaseOnEntities(DataForcesAndSourcesCore &data) const
Get field approximation space and base on entities.
Definition: ForcesAndSourcesCore.cpp:971
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:178
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:123
MoFEM::Types::VectorInt
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:73
MoFEM::ForcesAndSourcesCore::getMaxColOrder
int getMaxColOrder() const
Get max order of approximation for field in columns.
Definition: ForcesAndSourcesCore.cpp:142
MoFEM::DataForcesAndSourcesCore::spacesOnEntities
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
Definition: DataStructures.hpp:1034
MoFEM::VectorDofs
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
Definition: DataStructures.hpp:35
MoFEM::DataForcesAndSourcesCore::bAse
std::bitset< LASTBASE > bAse
bases on element
Definition: DataStructures.hpp:1030
MoFEM::Types::BitFieldId
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition: Types.hpp:53
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:150
cblas_ddot
double cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY)
Definition: cblas_ddot.c:12
MoFEM::ForcesAndSourcesCore::createDataOnElement
MoFEMErrorCode createDataOnElement()
Create a entity data on element object.
Definition: ForcesAndSourcesCore.cpp:1256
MoFEM::ContactPrismElementForcesAndSourcesCore::UserDataOperator::FACESLAVESLAVE
@ FACESLAVESLAVE
Definition: ContactPrismElementForcesAndSourcesCore.hpp:94
MoFEM::ContactPrismElementForcesAndSourcesCore::tangentMasterTwo
VectorDouble tangentMasterTwo
Definition: ContactPrismElementForcesAndSourcesCore.hpp:271
MoFEM::ForcesAndSourcesCore::opPtrVector
boost::ptr_vector< UserDataOperator > opPtrVector
Vector of finite element users data operators.
Definition: ForcesAndSourcesCore.hpp:816
MoFEM::ContactPrismElementForcesAndSourcesCore::getNodesIndices
MoFEMErrorCode getNodesIndices(const std::string field_name, FieldEntity_vector_view &ents_field, VectorInt &master_nodes_indices, VectorInt &master_local_nodes_indices, VectorInt &slave_nodes_indices, VectorInt &slave_local_nodes_indices, EXTRACTOR &&extractor) const
function that gets nodes indices.
Definition: ContactPrismElementForcesAndSourcesCore.cpp:1120
MoFEM::ContactPrismElementForcesAndSourcesCore::dataHdivSlave
DataForcesAndSourcesCore & dataHdivSlave
Definition: ContactPrismElementForcesAndSourcesCore.hpp:305
i
FTensor::Index< 'i', 3 > i
Definition: matrix_function.cpp:18
MoFEM::FieldEntity_vector_view
std::vector< boost::weak_ptr< FieldEntity > > FieldEntity_vector_view
Definition: FieldEntsMultiIndices.hpp:449
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEM::FEMethod::getRowFieldEnts
auto & getRowFieldEnts() const
Definition: LoopMethods.hpp:411
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:179
MoFEM::Tools::getTriNormal
static MoFEMErrorCode getTriNormal(const double *coords, double *normal)
Get the Tri Normal objectGet triangle normal.
Definition: Tools.cpp:258
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:124
MoFEM::Types::VectorDouble
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
MoFEM::ContactPrismElementForcesAndSourcesCore::setDefaultGaussPts
MoFEMErrorCode setDefaultGaussPts(const int rule)
Definition: ContactPrismElementForcesAndSourcesCore.cpp:121
MoFEM::EntPolynomialBaseCtx
Class used to pass element data to calculate base functions on tet,triangle,edge.
Definition: EntPolynomialBaseCtx.hpp:43
MoFEM::ContactPrismElementForcesAndSourcesCore::getEntityIndices
MoFEMErrorCode getEntityIndices(DataForcesAndSourcesCore &master_data, DataForcesAndSourcesCore &slave_data, const std::string &field_name, FieldEntity_vector_view &ents_field, const EntityType type_lo, const EntityType type_hi, EXTRACTOR &&extractor) const
function that gets entity indices.
Definition: ContactPrismElementForcesAndSourcesCore.cpp:1026
cblas_dcopy
void cblas_dcopy(const int N, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_dcopy.c:11
MoFEM::ContactPrismElementForcesAndSourcesCore::UserDataOperator
default operator for Contact Prism element
Definition: ContactPrismElementForcesAndSourcesCore.hpp:46
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
ReactionDiffusionEquation::r
const double r
rate factor
Definition: reaction_diffusion.cpp:33
MoFEM::FEMethod::feName
std::string feName
Name of finite element.
Definition: LoopMethods.hpp:386
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
NOFIELD
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:176
cache
BlockParamData * cache
Definition: multifield_plasticity.cpp:81
Generated on Tue Jan 19 2021 13:52:53 for MoFEM by Doxygen 1.8.20 and hosted at University of Glasgow