v0.14.0
SurfacePressureComplexForLazy.cpp
Go to the documentation of this file.
1 /* \file SurfacePressureComplexForLazy.cpp
2  */
3 
4 
5 
6 #include <MoFEM.hpp>
7 
8 using namespace MoFEM;
9 #include <MethodForForceScaling.hpp>
10 #include <SurfacePressure.hpp>
11 #include <SurfacePressureComplexForLazy.hpp>
12 
13 extern "C" {
14 
15 MoFEMErrorCode Traction_hierarchical(int order, int *order_edge, double *N,
16  double *N_face, double *N_edge[],
17  double *t, double *t_edge[],
18  double *t_face, double *traction, int gg);
19 MoFEMErrorCode Fext_h_hierarchical(
20  int order, int *order_edge, double *N, double *N_face, double *N_edge[],
21  double *diffN, double *diffN_face, double *diffN_edge[], double *t,
22  double *t_edge[], double *t_face, double *dofs_x, double *dofs_x_edge[],
23  double *dofs_x_face, double *idofs_x, double *idofs_x_edge[],
24  double *idofs_x_face, double *Fext, double *Fext_egde[], double *Fext_face,
25  double *iFext, double *iFext_egde[], double *iFext_face, int g_dim,
26  const double *g_w);
27 MoFEMErrorCode KExt_hh_hierarchical(
28  double eps, int order, int *order_edge, double *N, double *N_face,
29  double *N_edge[], double *diffN, double *diffN_face, double *diffN_edge[],
30  double *t, double *t_edge[], double *t_face, double *dofs_x,
31  double *dofs_x_edge[], double *dofs_x_face, double *KExt_hh,
32  double *KExt_egdeh[3], double *KExt_faceh, int g_dim, const double *g_w);
33 MoFEMErrorCode KExt_hh_hierarchical_edge(
34  double eps, int order, int *order_edge, double *N, double *N_face,
35  double *N_edge[], double *diffN, double *diffN_face, double *diffN_edge[],
36  double *t, double *t_edge[], double *t_face, double *dofs_x,
37  double *dofs_x_edge[], double *dofs_x_face, double *Khext_edge[3],
38  double *KExt_edgeegde[3][3], double *KExt_faceedge[3], int g_dim,
39  const double *g_w);
40 MoFEMErrorCode
41 KExt_hh_hierarchical_face(double eps, int order, int *order_edge, double *N,
42  double *N_face, double *N_edge[], double *diffN,
43  double *diffN_face, double *diffN_edge[], double *t,
44  double *t_edge[], double *t_face, double *dofs_x,
45  double *dofs_x_edge[], double *dofs_x_face,
46  double *KExt_hface, double *KExt_egdeface[3],
47  double *KExt_faceface, int g_dim, const double *g_w);
48 
49 void tetcircumcenter_tp(double a[3], double b[3], double c[3], double d[3],
50  double circumcenter[3], double *xi, double *eta,
51  double *zeta);
52 void tricircumcenter3d_tp(double a[3], double b[3], double c[3],
53  double circumcenter[3], double *xi, double *eta);
54 }
55 
56 NeumannForcesSurfaceComplexForLazy::AuxMethodSpatial::AuxMethodSpatial(
57  const string &field_name, MyTriangleSpatialFE *my_ptr, const char type)
58  : FaceElementForcesAndSourcesCore::UserDataOperator(field_name, type),
59  myPtr(my_ptr) {}
60 
61 NeumannForcesSurfaceComplexForLazy::AuxMethodMaterial::AuxMethodMaterial(
62  const string &field_name, MyTriangleSpatialFE *my_ptr, const char type)
63  : FaceElementForcesAndSourcesCore::UserDataOperator(field_name, type),
64  myPtr(my_ptr){};
65 
66 MoFEMErrorCode NeumannForcesSurfaceComplexForLazy::AuxMethodSpatial::doWork(
67  int side, EntityType type, EntitiesFieldData::EntData &data) {
68  MoFEMFunctionBegin;
69 
70  switch (type) {
71  case MBVERTEX: {
72  if (data.getFieldData().size() != 9) {
73  SETERRQ2(PETSC_COMM_SELF, 1,
74  "it should be 9 dofs on vertices but is %d of field < %s >",
75  data.getFieldData().size(), rowFieldName.c_str());
76  }
77  myPtr->N = &*data.getN().data().begin();
78  myPtr->diffN = &*data.getDiffN().data().begin();
79  myPtr->dOfs_x.resize(data.getFieldData().size());
80  ublas::noalias(myPtr->dOfs_x) = data.getFieldData();
81  myPtr->dofs_x = &*myPtr->dOfs_x.data().begin();
82  myPtr->dOfs_x_indices.resize(data.getIndices().size());
83  ublas::noalias(myPtr->dOfs_x_indices) = data.getIndices();
84  myPtr->dofs_x_indices = &*myPtr->dOfs_x_indices.data().begin();
85  } break;
86  case MBEDGE: {
87  myPtr->order_edge[side] = data.getOrder();
88  myPtr->N_edge[side] = &*data.getN().data().begin();
89  myPtr->diffN_edge[side] = &*data.getDiffN().data().begin();
90  myPtr->dOfs_x_edge.resize(3);
91  myPtr->dOfs_x_edge[side].resize(data.getFieldData().size());
92  myPtr->dofs_x_edge[side] = &*myPtr->dOfs_x_edge[side].data().begin();
93  myPtr->dOfs_x_edge_indices.resize(3);
94  myPtr->dOfs_x_edge_indices[side].resize(data.getIndices().size());
95  ublas::noalias(myPtr->dOfs_x_edge_indices[side]) = data.getIndices();
96  myPtr->dofs_x_edge_indices[side] =
97  &*myPtr->dOfs_x_edge_indices[side].data().begin();
98  } break;
99  case MBTRI: {
100  myPtr->order_face = data.getOrder();
101  myPtr->N_face = &*data.getN().data().begin();
102  myPtr->diffN_face = &*data.getDiffN().data().begin();
103  myPtr->dOfs_x_face.resize(data.getFieldData().size());
104  ublas::noalias(myPtr->dOfs_x_face) = data.getFieldData();
105  myPtr->dofs_x_face = &*myPtr->dOfs_x_face.data().begin();
106  myPtr->dOfs_x_face_indices.resize(data.getIndices().size());
107  ublas::noalias(myPtr->dOfs_x_face_indices) = data.getIndices();
108  myPtr->dofs_x_face_indices = &*myPtr->dOfs_x_face_indices.data().begin();
109  } break;
110  default:
111  SETERRQ(PETSC_COMM_SELF, 1, "unknown entity type");
112  }
113 
114  MoFEMFunctionReturn(0);
115 }
116 
117 MoFEMErrorCode NeumannForcesSurfaceComplexForLazy::AuxMethodMaterial::doWork(
118  int side, EntityType type, EntitiesFieldData::EntData &data) {
119  MoFEMFunctionBegin;
120 
121  switch (type) {
122  case MBVERTEX: {
123  if (data.getFieldData().size() != 9) {
124  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
125  "it should be 9 dofs on vertices but is %d",
126  data.getFieldData().size());
127  }
128  if (data.getN().size2() != 3) {
129  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
130  "it should 3 shape functions for 3 nodes");
131  }
132  myPtr->N = &*data.getN().data().begin();
133  myPtr->diffN = &*data.getDiffN().data().begin();
134  myPtr->dOfs_X_indices.resize(data.getIndices().size());
135  ublas::noalias(myPtr->dOfs_X_indices) = data.getIndices();
136  myPtr->dofs_X_indices = &*myPtr->dOfs_X_indices.data().begin();
137  myPtr->dOfs_X.resize(data.getFieldData().size());
138  ublas::noalias(myPtr->dOfs_X) = data.getFieldData();
139  myPtr->dofs_X = &*myPtr->dOfs_X.data().begin();
140  } break;
141  case MBEDGE: {
142  myPtr->order_edge_material[side] = data.getOrder();
143  myPtr->dOfs_X_edge.resize(3);
144  myPtr->dOfs_X_edge[side].resize(data.getFieldData().size());
145  ublas::noalias(myPtr->dOfs_X_edge[side]) = data.getFieldData();
146  myPtr->dofs_X_edge[side] = &*myPtr->dOfs_X_edge[side].data().begin();
147  } break;
148  case MBTRI: {
149  myPtr->order_face_material = data.getOrder();
150  myPtr->dOfs_X_face.resize(data.getFieldData().size());
151  ublas::noalias(myPtr->dOfs_X_face) = data.getFieldData();
152  myPtr->dofs_X_face = &*myPtr->dOfs_X_face.data().begin();
153  } break;
154  default:
155  SETERRQ(PETSC_COMM_SELF, 1, "unknown entity type");
156  }
157 
158  MoFEMFunctionReturn(0);
159 }
160 
161 NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::MyTriangleSpatialFE(
162  MoFEM::Interface &_mField, Mat _Aij, Vec &_F, double *scale_lhs,
163  double *scale_rhs, std::string spatial_field_name,
164  std::string mat_field_name)
165  : FaceElementForcesAndSourcesCore(_mField), sCaleLhs(scale_lhs),
166  sCaleRhs(scale_rhs), typeOfForces(CONSERVATIVE), eps(1e-8), uSeF(false) {
167 
168  meshPositionsFieldName = "NoNE";
169  methodsOp.clear();
170 
171  Aij = _Aij;
172  F = _F;
173 
174  snes_B = _Aij;
175  snes_f = _F;
176 
177  if (mField.check_field(mat_field_name)) {
178  getOpPtrVector().push_back(new AuxMethodMaterial(
179  mat_field_name, this, ForcesAndSourcesCore::UserDataOperator::OPROW));
180  }
181  getOpPtrVector().push_back(new AuxMethodSpatial(
182  spatial_field_name, this, ForcesAndSourcesCore::UserDataOperator::OPROW));
183 }
184 
185 MoFEMErrorCode NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::rHs() {
186  MoFEMFunctionBegin;
187 
188  auto &dataH1 = *dataOnElement[H1];
189 
190  fExtNode.resize(9);
191  fExtFace.resize(dataH1.dataOnEntities[MBTRI][0].getFieldData().size());
192  fExtEdge.resize(3);
193  for (int ee = 0; ee < 3; ee++) {
194  int nb_edge_dofs = dOfs_x_edge_indices[ee].size();
195  if (nb_edge_dofs > 0) {
196  fExtEdge[ee].resize(nb_edge_dofs);
197  Fext_edge[ee] = &*fExtEdge[ee].data().begin();
198  } else {
199  Fext_edge[ee] = NULL;
200  }
201  }
202 
203  switch (typeOfForces) {
204  case CONSERVATIVE:
205  CHKERR Fext_h_hierarchical(
206  order_face, order_edge, // 2
207  N, N_face, N_edge, diffN, diffN_face, diffN_edge, // 8
208  t_loc, NULL, NULL, // 11
209  dofs_x, dofs_x_edge, dofs_x_face, // 14
210  NULL, NULL, NULL, // 17
211  &*fExtNode.data().begin(), Fext_edge, &*fExtFace.data().begin(), // 20
212  NULL, NULL, NULL, // 23
213  gaussPts.size2(), &gaussPts(2, 0));
214  break;
215  case NONCONSERVATIVE:
216  for (int ee = 0; ee < 3; ee++) {
217  dOfs_X_edge.resize(3);
218  unsigned int s = dOfs_X_edge[ee].size();
219  dOfs_X_edge[ee].resize(dOfs_x_edge[ee].size(), true);
220  for (; s < dOfs_X_edge[ee].size(); s++) {
221  dOfs_X_edge[ee][s] = 0;
222  }
223  dofs_X_edge[ee] = &*dOfs_X_edge[ee].data().begin();
224  }
225  unsigned int s = dOfs_X_face.size();
226  dOfs_X_face.resize(dOfs_x_face.size(), true);
227  for (; s < dOfs_X_face.size(); s++) {
228  dOfs_X_face[s] = 0;
229  }
230  dofs_X_face = &*dOfs_X_face.data().begin();
231 
232  CHKERR Fext_h_hierarchical(
233  order_face, order_edge, // 2
234  N, N_face, N_edge, diffN, diffN_face, diffN_edge, // 8
235  t_loc, NULL, NULL, // 11
236  dofs_X, dofs_X_edge, dofs_X_face, // 14
237  NULL, NULL, NULL, // 17
238  &*fExtNode.data().begin(), Fext_edge, &*fExtFace.data().begin(), // 20
239  NULL, NULL, NULL, // 23
240  gaussPts.size2(), &gaussPts(2, 0));
241  break;
242  }
243 
244  Vec f = snes_f;
245  if (uSeF)
246  f = F;
247 
248  CHKERR VecSetValues(f, 9, dofs_x_indices, &*fExtNode.data().begin(),
249  ADD_VALUES);
250  if (dOfs_x_face_indices.size() > 0) {
251  CHKERR VecSetValues(f, dOfs_x_face_indices.size(), dofs_x_face_indices,
252  &*fExtFace.data().begin(), ADD_VALUES);
253  }
254  for (int ee = 0; ee < 3; ee++) {
255  if (dOfs_x_edge_indices[ee].size() > 0) {
256  CHKERR VecSetValues(f, dOfs_x_edge_indices[ee].size(),
257  dofs_x_edge_indices[ee], Fext_edge[ee], ADD_VALUES);
258  }
259  }
260 
261  MoFEMFunctionReturn(0);
262 }
263 
264 MoFEMErrorCode NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::lHs() {
265  MoFEMFunctionBegin;
266 
267  if (typeOfForces == NONCONSERVATIVE) {
268  MoFEMFunctionReturnHot(0);
269  }
270 
271  auto &dataH1 = *dataOnElement[H1];
272 
273  double center[3];
274  tricircumcenter3d_tp(&coords.data()[0], &coords.data()[3], &coords.data()[6],
275  center, NULL, NULL);
276  cblas_daxpy(3, -1, &coords.data()[0], 1, center, 1);
277  double r = cblas_dnrm2(3, center, 1);
278 
279  kExtNodeNode.resize(9, 9);
280  kExtEdgeNode.resize(3);
281  for (int ee = 0; ee < 3; ee++) {
282  kExtEdgeNode[ee].resize(dOfs_x_edge_indices[ee].size(), 9);
283  Kext_edge_node[ee] = &*kExtEdgeNode[ee].data().begin();
284  }
285  kExtFaceNode.resize(dOfs_x_face_indices.size(), 9);
286  CHKERR KExt_hh_hierarchical(
287  r * eps, order_face, order_edge, N, N_face, N_edge, diffN, diffN_face,
288  diffN_edge, t_loc, NULL, NULL, dofs_x, dofs_x_edge, dofs_x_face,
289  &*kExtNodeNode.data().begin(), Kext_edge_node,
290  &*kExtFaceNode.data().begin(), gaussPts.size2(), &gaussPts(2, 0));
291  // cerr << kExtNodeNode << endl;
292  CHKERR MatSetValues(snes_B, 9, dofs_x_indices, 9, dofs_x_indices,
293  &*kExtNodeNode.data().begin(), ADD_VALUES);
294  CHKERR MatSetValues(snes_B, kExtFaceNode.size1(), dofs_x_face_indices, 9,
295  dofs_x_indices, &*kExtFaceNode.data().begin(),
296  ADD_VALUES);
297  // cerr << kExtFaceNode << endl;
298  for (int ee = 0; ee < 3; ee++) {
299  // cerr << kExtEdgeNode[ee] << endl;
300  CHKERR MatSetValues(snes_B, kExtEdgeNode[ee].size1(),
301  dofs_x_edge_indices[ee], 9, dofs_x_indices,
302  Kext_edge_node[ee], ADD_VALUES);
303  }
304 
305  kExtNodeFace.resize(9, dOfs_x_face_indices.size());
306  kExtEdgeFace.resize(3);
307  for (int ee = 0; ee < 3; ee++) {
308  kExtEdgeFace[ee].resize(
309  dOfs_x_edge_indices[ee].size(),
310  dataH1.dataOnEntities[MBTRI][0].getIndices().size());
311  Kext_edge_face[ee] = &*kExtEdgeFace[ee].data().begin();
312  }
313  kExtFaceFace.resize(dOfs_x_face_indices.size(), dOfs_x_face_indices.size());
314  CHKERR KExt_hh_hierarchical_face(
315  r * eps, order_face, order_edge, N, N_face, N_edge, diffN, diffN_face,
316  diffN_edge, t_loc, NULL, NULL, dofs_x, dofs_x_edge, dofs_x_face,
317  &*kExtNodeFace.data().begin(), Kext_edge_face,
318  &*kExtFaceFace.data().begin(), gaussPts.size2(), &gaussPts(2, 0));
319  // cerr << "kExtNodeFace " << kExtNodeFace << endl;
320  // cerr << "kExtFaceFace " << kExtFaceFace << endl;
321  CHKERR MatSetValues(snes_B, 9, dofs_x_indices, kExtNodeFace.size2(),
322  dofs_x_face_indices, &*kExtNodeFace.data().begin(),
323  ADD_VALUES);
324  CHKERR MatSetValues(snes_B, kExtFaceFace.size1(), dofs_x_face_indices,
325  kExtFaceFace.size2(), dofs_x_face_indices,
326  &*kExtFaceFace.data().begin(), ADD_VALUES);
327  for (int ee = 0; ee < 3; ee++) {
328  // cerr << "kExtEdgeFace " << kExtEdgeFace[ee] << endl;
329  CHKERR MatSetValues(snes_B, kExtEdgeFace[ee].size1(),
330  dofs_x_edge_indices[ee], kExtFaceFace.size2(),
331  dofs_x_face_indices, Kext_edge_face[ee], ADD_VALUES);
332  }
333 
334  kExtFaceEdge.resize(3);
335  kExtNodeEdge.resize(3);
336  kExtEdgeEdge.resize(3, 3);
337  for (int ee = 0; ee < 3; ee++) {
338  if (dOfs_x_edge_indices[ee].size() !=
339  (unsigned int)(3 * NBEDGE_H1(order_edge[ee]))) {
340  SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
341  }
342  kExtFaceEdge[ee].resize(dOfs_x_face_indices.size(),
343  dOfs_x_edge_indices[ee].size());
344  kExtNodeEdge[ee].resize(9, dOfs_x_edge_indices[ee].size());
345  Kext_node_edge[ee] = &*kExtNodeEdge[ee].data().begin();
346  Kext_face_edge[ee] = &*kExtFaceEdge[ee].data().begin();
347  for (int EE = 0; EE < 3; EE++) {
348  kExtEdgeEdge(EE, ee).resize(dOfs_x_edge_indices[EE].size(),
349  dOfs_x_edge_indices[ee].size());
350  Kext_edge_edge[EE][ee] = &*kExtEdgeEdge(EE, ee).data().begin();
351  }
352  }
353  CHKERR KExt_hh_hierarchical_edge(
354  r * eps, order_face, order_edge, N, N_face, N_edge, diffN, diffN_face,
355  diffN_edge, t_loc, NULL, NULL, dofs_x, dofs_x_edge, dofs_x_face,
356  Kext_node_edge, Kext_edge_edge, Kext_face_edge, gaussPts.size2(),
357  &gaussPts(2, 0));
358  for (int ee = 0; ee < 3; ee++) {
359  CHKERR MatSetValues(snes_B, kExtFaceEdge[ee].size1(), dofs_x_face_indices,
360  kExtFaceEdge[ee].size2(), dofs_x_edge_indices[ee],
361  &*kExtFaceEdge[ee].data().begin(), ADD_VALUES);
362  CHKERR MatSetValues(snes_B, 9, dofs_x_indices, kExtNodeEdge[ee].size2(),
363  dofs_x_edge_indices[ee],
364  &*kExtNodeEdge[ee].data().begin(), ADD_VALUES);
365  for (int EE = 0; EE < 3; EE++) {
366  CHKERR MatSetValues(snes_B, kExtEdgeEdge(EE, ee).size1(),
367  dofs_x_edge_indices[EE], kExtEdgeEdge(EE, ee).size2(),
368  dofs_x_edge_indices[ee], Kext_edge_edge[EE][ee],
369  ADD_VALUES);
370  }
371  }
372 
373  MoFEMFunctionReturn(0);
374 }
375 
376 MoFEMErrorCode NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::
377  reBaseToFaceLoocalCoordSystem(MatrixDouble &t_glob_nodal) {
378  MoFEMFunctionBeginHot;
379  double s1[3], s2[3], normal[3], q[9];
380  CHKERR ShapeFaceBaseMBTRI(diffN, &*coords.data().begin(), normal, s1, s2);
381  double nrm2_normal = cblas_dnrm2(3, normal, 1);
382  cblas_dscal(3, 1. / nrm2_normal, normal, 1);
383  cblas_dcopy(3, s1, 1, &q[0], 1);
384  cblas_dcopy(3, s2, 1, &q[3], 1);
385  cblas_dcopy(3, normal, 1, &q[6], 1);
386  __CLPK_integer info;
387  __CLPK_integer ipiv[3];
388  info = lapack_dgesv(3, 3, q, 3, ipiv, &*t_glob_nodal.data().begin(), 3);
389  if (info != 0) {
390  SETERRQ1(PETSC_COMM_SELF, 1, "error solve dgesv info = %d", info);
391  }
392  MoFEMFunctionReturnHot(0);
393 }
394 
395 MoFEMErrorCode
396 NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::calcTraction() {
397  MoFEMFunctionBeginHot;
398  EntityHandle ent = numeredEntFiniteElementPtr->getEnt();
399  map<int, bCPressure>::iterator mip = mapPressure.begin();
400  tLoc.resize(3);
401  tLoc[0] = tLoc[1] = tLoc[2] = 0;
402  for (; mip != mapPressure.end(); mip++) {
403  if (mip->second.tRis.find(ent) != mip->second.tRis.end()) {
404  tLoc[2] -= mip->second.data.data.value1;
405  }
406  }
407  tLocNodal.resize(3, 3);
408  for (int nn = 0; nn < 3; nn++) {
409  for (int dd = 0; dd < 3; dd++) {
410  tLocNodal(nn, dd) = tLoc[dd];
411  }
412  }
413 
414  map<int, bCForce>::iterator mif = mapForce.begin();
415  for (; mif != mapForce.end(); mif++) {
416  if (mif->second.tRis.find(ent) != mif->second.tRis.end()) {
417  tGlob.resize(3);
418  tGlob[0] = mif->second.data.data.value3;
419  tGlob[1] = mif->second.data.data.value4;
420  tGlob[2] = mif->second.data.data.value5;
421  tGlob *= mif->second.data.data.value1;
422  tGlobNodal.resize(3, 3);
423  for (int nn = 0; nn < 3; nn++) {
424  for (int dd = 0; dd < 3; dd++) {
425  tGlobNodal(nn, dd) = tGlob[dd];
426  }
427  }
428  CHKERR reBaseToFaceLoocalCoordSystem(tGlobNodal);
429  tLocNodal += tGlobNodal;
430  }
431  }
432 
433  VectorDouble scale(1, 1);
434  CHKERR MethodForForceScaling::applyScale(this, methodsOp, scale);
435  tLocNodal *= scale[0];
436 
437  t_loc = &*tLocNodal.data().begin();
438 
439  MoFEMFunctionReturnHot(0);
440 }
441 
442 MoFEMErrorCode
443 NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::preProcess() {
444  MoFEMFunctionBeginHot;
445 
446  CHKERR PetscOptionsBegin(mField.get_comm(), "",
447  "Surface Pressure (complex for lazy)", "none");
448  PetscBool is_conservative = PETSC_TRUE;
449  CHKERR PetscOptionsBool("-is_conservative_force", "is conservative force", "",
450  PETSC_TRUE, &is_conservative, PETSC_NULL);
451  if (is_conservative == PETSC_FALSE) {
452  typeOfForces = NONCONSERVATIVE;
453  }
454  ierr = PetscOptionsEnd();
455  CHKERRG(ierr);
456 
457  MoFEMFunctionReturnHot(0);
458 }
459 
460 MoFEMErrorCode NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::
461 operator()() {
462  MoFEMFunctionBeginHot;
463 
464  {
465 
466  {
467 
468  dofs_X = &*coords.data().begin();
469  for (int ee = 0; ee < 3; ee++) {
470  dofs_X_edge[ee] = NULL;
471  idofs_X_edge[ee] = NULL;
472  order_edge_material[ee] = 0;
473  }
474  dofs_X_face = NULL;
475  idofs_X_face = NULL;
476  order_face_material = 0;
477 
478  dofs_x = &*coords.data().begin();
479  idofs_x = NULL;
480  for (int ee = 0; ee < 3; ee++) {
481  order_edge[ee] = 0;
482  N_edge[ee] = NULL;
483  diffN_edge[ee] = NULL;
484  dofs_x_edge[ee] = NULL;
485  idofs_x_edge[ee] = NULL;
486  }
487  order_face = 0;
488  N_face = NULL;
489  diffN_face = NULL;
490  dofs_x_face = NULL;
491  idofs_x_face = NULL;
492  }
493 
494  CHKERR FaceElementForcesAndSourcesCore::operator()();
495  CHKERR calcTraction();
496 
497  switch (snes_ctx) {
498  case CTX_SNESNONE:
499  case CTX_SNESSETFUNCTION: {
500  tLocNodal *= *sCaleRhs;
501  // cerr << "sCaleRhs " << *sCaleRhs << endl;
502  // cerr << tLocNodal << endl;
503  CHKERR rHs();
504  } break;
505  case CTX_SNESSETJACOBIAN: {
506  tLocNodal *= *sCaleLhs;
507  CHKERR lHs();
508  } break;
509  }
510  }
511 
512  MoFEMFunctionReturnHot(0);
513 }
514 
515 MoFEMErrorCode
516 NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::addForce(int ms_id) {
517  MeshsetsManager *mesh_manager_ptr;
518  const CubitMeshSets *cubit_meshset_ptr;
519 
520  MoFEMFunctionBeginHot;
521  CHKERR mField.getInterface(mesh_manager_ptr);
522  CHKERR mesh_manager_ptr->getCubitMeshsetPtr(ms_id, NODESET,
523  &cubit_meshset_ptr);
524  CHKERR cubit_meshset_ptr->getBcDataStructure(mapForce[ms_id].data);
525  CHKERR mField.get_moab().get_entities_by_type(
526  cubit_meshset_ptr->meshset, MBTRI, mapForce[ms_id].tRis, true);
527  MoFEMFunctionReturnHot(0);
528 }
529 
530 MoFEMErrorCode
531 NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::addPressure(
532  int ms_id) {
533  MeshsetsManager *mesh_manager_ptr;
534  const CubitMeshSets *cubit_meshset_ptr;
535 
536  MoFEMFunctionBeginHot;
537  CHKERR mField.getInterface(mesh_manager_ptr);
538  CHKERR mesh_manager_ptr->getCubitMeshsetPtr(ms_id, SIDESET,
539  &cubit_meshset_ptr);
540  CHKERR cubit_meshset_ptr->getBcDataStructure(mapPressure[ms_id].data);
541  CHKERR mField.get_moab().get_entities_by_type(
542  cubit_meshset_ptr->meshset, MBTRI, mapPressure[ms_id].tRis, true);
543  MoFEMFunctionReturnHot(0);
544 }
545 
546 NeumannForcesSurfaceComplexForLazy::NeumannForcesSurfaceComplexForLazy(
547  MoFEM::Interface &m_field, Mat _Aij, Vec _F, std::string spatial_field_name,
548  std::string material_field_name)
549  : mField(m_field), feSpatial(m_field, _Aij, _F, NULL, NULL,
550  spatial_field_name, material_field_name),
551  spatialField(spatial_field_name), materialField(material_field_name) {
552 
553  double def_scale = 1.;
554  const EntityHandle root_meshset = mField.get_moab().get_root_set();
555  CHKERR mField.get_moab().tag_get_handle(
556  "_LoadFactor_Scale_", 1, MB_TYPE_DOUBLE, thScale,
557  MB_TAG_CREAT | MB_TAG_EXCL | MB_TAG_MESH, &def_scale);
558  if (rval == MB_ALREADY_ALLOCATED) {
559  rval = mField.get_moab().tag_get_by_ptr(thScale, &root_meshset, 1,
560  (const void **)&sCale);
561  MOAB_THROW(rval);
562  } else {
563  MOAB_THROW(rval);
564  rval =
565  mField.get_moab().tag_set_data(thScale, &root_meshset, 1, &def_scale);
566  MOAB_THROW(rval);
567  rval = mField.get_moab().tag_get_by_ptr(thScale, &root_meshset, 1,
568  (const void **)&sCale);
569  MOAB_THROW(rval);
570  }
571 
572  feSpatial.sCaleLhs = sCale;
573  feSpatial.sCaleRhs = sCale;
574 }
575 
576 NeumannForcesSurfaceComplexForLazy::NeumannForcesSurfaceComplexForLazy(
577  MoFEM::Interface &m_field, Mat _Aij, Vec _F, double *scale_lhs,
578  double *scale_rhs, std::string spatial_field_name,
579  std::string material_field_name)
580  : mField(m_field), feSpatial(m_field, _Aij, _F, scale_lhs, scale_rhs,
581  spatial_field_name, material_field_name),
582  spatialField(spatial_field_name), materialField(material_field_name) {}
583 
584 extern "C" {
585 
586 // External forces
587 MoFEMErrorCode Traction_hierarchical(int order, int *order_edge, double *N,
588  double *N_face, double *N_edge[],
589  double *t, double *t_edge[],
590  double *t_face, double *traction, int gg) {
591  MoFEMFunctionBeginHot;
592  int dd, ee;
593  for (dd = 0; dd < 3; dd++)
594  traction[dd] = cblas_ddot(3, &N[gg * 3], 1, &t[dd], 3);
595  if (t_face != NULL) {
596  int nb_dofs_face = NBFACETRI_H1(order);
597  if (nb_dofs_face > 0) {
598  for (dd = 0; dd < 3; dd++)
599  traction[dd] += cblas_ddot(nb_dofs_face, &N_face[gg * nb_dofs_face], 1,
600  &t_face[dd], 3);
601  }
602  }
603  if (t_edge != NULL) {
604  ee = 0;
605  for (; ee < 3; ee++) {
606  if (t_edge[ee] == NULL)
607  continue;
608  int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
609  if (nb_dofs_edge > 0) {
610  for (dd = 0; dd < 3; dd++) {
611  traction[dd] +=
612  cblas_ddot(nb_dofs_edge, &(N_edge[ee][gg * nb_dofs_edge]), 1,
613  &(t_edge[ee][dd]), 3);
614  }
615  }
616  }
617  }
618  MoFEMFunctionReturnHot(0);
619 }
620 MoFEMErrorCode Fext_h_hierarchical(
621  int order, int *order_edge, double *N, double *N_face, double *N_edge[],
622  double *diffN, double *diffN_face, double *diffN_edge[], double *t,
623  double *t_edge[], double *t_face, double *dofs_x, double *dofs_x_edge[],
624  double *dofs_x_face, double *idofs_x, double *idofs_x_edge[],
625  double *idofs_x_face, double *Fext, double *Fext_edge[], double *Fext_face,
626  double *iFext, double *iFext_edge[], double *iFext_face, int g_dim,
627  const double *g_w) {
628  MoFEMFunctionBeginHot;
629  int dd, nn, ee, gg;
630  if (Fext != NULL)
631  bzero(Fext, 9 * sizeof(double));
632  if (iFext != NULL)
633  bzero(iFext, 9 * sizeof(double));
634  ee = 0;
635  for (; ee < 3; ee++) {
636  int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
637  if (nb_dofs_edge == 0)
638  continue;
639  if (Fext_edge != NULL)
640  bzero(Fext_edge[ee], 3 * nb_dofs_edge * sizeof(double));
641  if (iFext_edge != NULL)
642  bzero(iFext_edge[ee], 3 * nb_dofs_edge * sizeof(double));
643  }
644  int nb_dofs_face = NBFACETRI_H1(order);
645  if (nb_dofs_face != 0) {
646  if (Fext_face != NULL)
647  bzero(Fext_face, 3 * nb_dofs_face * sizeof(double));
648  if (iFext_face != NULL)
649  bzero(iFext_face, 3 * nb_dofs_face * sizeof(double));
650  }
651  gg = 0;
652  for (; gg < g_dim; gg++) {
653  double traction[3] = {0, 0, 0};
654  CHKERR Traction_hierarchical(order, order_edge, N, N_face, N_edge, t,
655  t_edge, t_face, traction, gg);
656  __CLPK_doublecomplex xnormal[3], xs1[3], xs2[3];
657  CHKERR Normal_hierarchical(
658  // FIXME: order of tractions approximation could be different for field
659  // approx.
660  order, order_edge, order, order_edge, diffN, diffN_face, diffN_edge,
661  dofs_x, dofs_x_edge, dofs_x_face, idofs_x, idofs_x_edge, idofs_x_face,
662  xnormal, xs1, xs2, gg);
663  CHKERR Base_scale(xnormal, xs1, xs2);
664  double normal_real[3], s1_real[3], s2_real[3];
665  double normal_imag[3], s1_imag[3], s2_imag[3];
666  for (dd = 0; dd < 3; dd++) {
667  normal_real[dd] = 0.5 * xnormal[dd].r;
668  normal_imag[dd] = 0.5 * xnormal[dd].i;
669  s1_real[dd] = 0.5 * xs1[dd].r;
670  s1_imag[dd] = 0.5 * xs1[dd].i;
671  s2_real[dd] = 0.5 * xs2[dd].r;
672  s2_imag[dd] = 0.5 * xs2[dd].i;
673  }
674  nn = 0;
675  for (; nn < 3; nn++) {
676  if (Fext != NULL)
677  for (dd = 0; dd < 3; dd++) {
678  // fprintf(stderr,"%d %f %f %f %f %f %f\n",
679  // gg,g_w[gg],N[3*gg+nn],normal_real[dd],
680  // traction[0],traction[1],traction[2]);
681  Fext[3 * nn + dd] +=
682  g_w[gg] * N[3 * gg + nn] * normal_real[dd] * traction[2];
683  Fext[3 * nn + dd] +=
684  g_w[gg] * N[3 * gg + nn] * s1_real[dd] * traction[0];
685  Fext[3 * nn + dd] +=
686  g_w[gg] * N[3 * gg + nn] * s2_real[dd] * traction[1];
687  }
688  if (iFext != NULL)
689  for (dd = 0; dd < 3; dd++) {
690  iFext[3 * nn + dd] +=
691  g_w[gg] * N[3 * gg + nn] * normal_imag[dd] * traction[2];
692  iFext[3 * nn + dd] +=
693  g_w[gg] * N[3 * gg + nn] * s1_imag[dd] * traction[0];
694  iFext[3 * nn + dd] +=
695  g_w[gg] * N[3 * gg + nn] * s2_imag[dd] * traction[1];
696  }
697  }
698  ee = 0;
699  for (; ee < 3; ee++) {
700  int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
701  if (nb_dofs_edge == 0)
702  continue;
703  int nn = 0;
704  for (; nn < nb_dofs_edge; nn++) {
705  if (Fext_edge != NULL)
706  for (dd = 0; dd < 3; dd++) {
707  Fext_edge[ee][3 * nn + dd] += g_w[gg] *
708  N_edge[ee][gg * nb_dofs_edge + nn] *
709  normal_real[dd] * traction[2];
710  Fext_edge[ee][3 * nn + dd] += g_w[gg] *
711  N_edge[ee][gg * nb_dofs_edge + nn] *
712  s1_real[dd] * traction[0];
713  Fext_edge[ee][3 * nn + dd] += g_w[gg] *
714  N_edge[ee][gg * nb_dofs_edge + nn] *
715  s2_real[dd] * traction[1];
716  }
717  if (iFext_edge != NULL) {
718  for (dd = 0; dd < 3; dd++) {
719  iFext_edge[ee][3 * nn + dd] += g_w[gg] *
720  N_edge[ee][gg * nb_dofs_edge + nn] *
721  normal_imag[dd] * traction[2];
722  iFext_edge[ee][3 * nn + dd] += g_w[gg] *
723  N_edge[ee][gg * nb_dofs_edge + nn] *
724  s1_imag[dd] * traction[0];
725  iFext_edge[ee][3 * nn + dd] += g_w[gg] *
726  N_edge[ee][gg * nb_dofs_edge + nn] *
727  s2_imag[dd] * traction[1];
728  }
729  }
730  }
731  }
732  if (nb_dofs_face != 0) {
733  nn = 0;
734  for (; nn < nb_dofs_face; nn++) {
735  if (Fext_face != NULL)
736  for (dd = 0; dd < 3; dd++) {
737  Fext_face[3 * nn + dd] += g_w[gg] * N_face[gg * nb_dofs_face + nn] *
738  normal_real[dd] * traction[2];
739  Fext_face[3 * nn + dd] += g_w[gg] * N_face[gg * nb_dofs_face + nn] *
740  s1_real[dd] * traction[0];
741  Fext_face[3 * nn + dd] += g_w[gg] * N_face[gg * nb_dofs_face + nn] *
742  s2_real[dd] * traction[1];
743  }
744  if (iFext_face != NULL)
745  for (dd = 0; dd < 3; dd++) {
746  iFext_face[3 * nn + dd] += g_w[gg] *
747  N_face[gg * nb_dofs_face + nn] *
748  normal_imag[dd] * traction[2];
749  iFext_face[3 * nn + dd] += g_w[gg] *
750  N_face[gg * nb_dofs_face + nn] *
751  s1_imag[dd] * traction[0];
752  iFext_face[3 * nn + dd] += g_w[gg] *
753  N_face[gg * nb_dofs_face + nn] *
754  s2_imag[dd] * traction[1];
755  }
756  }
757  }
758  }
759  MoFEMFunctionReturnHot(0);
760 }
761 MoFEMErrorCode KExt_hh_hierarchical(
762  double eps, int order, int *order_edge, double *N, double *N_face,
763  double *N_edge[], double *diffN, double *diffN_face, double *diffN_edge[],
764  double *t, double *t_edge[], double *t_face, double *dofs_x,
765  double *dofs_x_edge[], double *dofs_x_face, double *KExt_hh,
766  double *KExt_edgeh[], double *KExt_faceh, int g_dim, const double *g_w) {
767  MoFEMFunctionBeginHot;
768  int gg, dd, ii, nn, ee;
769  bzero(KExt_hh, 9 * 9 * sizeof(double));
770  if (KExt_edgeh != NULL) {
771  for (ee = 0; ee < 3; ee++) {
772  int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
773  bzero(KExt_edgeh[ee], 9 * 3 * nb_dofs_edge * sizeof(double));
774  }
775  }
776  int nb_dofs_face = NBFACETRI_H1(order);
777  if (KExt_faceh != NULL) {
778  bzero(KExt_faceh, 9 * 3 * nb_dofs_face * sizeof(double));
779  }
780  for (gg = 0; gg < g_dim; gg++) {
781  double traction[3] = {0, 0, 0};
782  CHKERR Traction_hierarchical(order, order_edge, N, N_face, N_edge, t,
783  t_edge, t_face, traction, gg);
784  //
785  __CLPK_doublecomplex xnormal[3], xs1[3], xs2[3];
786  double idofs_x[9];
787  for (ii = 0; ii < 9; ii++) {
788  bzero(idofs_x, 9 * sizeof(double));
789  idofs_x[ii] = eps;
790  CHKERR Normal_hierarchical(order, order_edge, // FIXME
791  order, order_edge, diffN, diffN_face,
792  diffN_edge, dofs_x, dofs_x_edge, dofs_x_face,
793  idofs_x, NULL, NULL, xnormal, xs1, xs2, gg);
794  CHKERR Base_scale(xnormal, xs1, xs2);
795  double normal_imag[3], s1_imag[3], s2_imag[3];
796  for (dd = 0; dd < 3; dd++) {
797  normal_imag[dd] = 0.5 * xnormal[dd].i / eps;
798  s1_imag[dd] = 0.5 * xs1[dd].i / eps;
799  s2_imag[dd] = 0.5 * xs2[dd].i / eps;
800  }
801  nn = 0;
802  for (; nn < 3; nn++) {
803  for (dd = 0; dd < 3; dd++) {
804  KExt_hh[ii + 9 * 3 * nn + 9 * dd] +=
805  g_w[gg] * N[3 * gg + nn] * normal_imag[dd] * traction[2];
806  KExt_hh[ii + 9 * 3 * nn + 9 * dd] +=
807  g_w[gg] * N[3 * gg + nn] * s1_imag[dd] * traction[0];
808  KExt_hh[ii + 9 * 3 * nn + 9 * dd] +=
809  g_w[gg] * N[3 * gg + nn] * s2_imag[dd] * traction[1];
810  }
811  }
812  if (KExt_edgeh != NULL) {
813  for (ee = 0; ee < 3; ee++) {
814  int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
815  for (nn = 0; nn < nb_dofs_edge; nn++) {
816  for (dd = 0; dd < 3; dd++) {
817  KExt_edgeh[ee][ii + 9 * 3 * nn + 9 * dd] +=
818  g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] *
819  normal_imag[dd] * traction[2];
820  KExt_edgeh[ee][ii + 9 * 3 * nn + 9 * dd] +=
821  g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] * s1_imag[dd] *
822  traction[0];
823  KExt_edgeh[ee][ii + 9 * 3 * nn + 9 * dd] +=
824  g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] * s2_imag[dd] *
825  traction[1];
826  }
827  }
828  }
829  }
830  if (KExt_faceh != NULL) {
831  for (nn = 0; nn < nb_dofs_face; nn++) {
832  for (dd = 0; dd < 3; dd++) {
833  KExt_faceh[ii + 3 * 9 * nn + 9 * dd] +=
834  g_w[gg] * N_face[nb_dofs_face * gg + nn] * normal_imag[dd] *
835  traction[2];
836  KExt_faceh[ii + 3 * 9 * nn + 9 * dd] +=
837  g_w[gg] * N_face[nb_dofs_face * gg + nn] * s1_imag[dd] *
838  traction[0];
839  KExt_faceh[ii + 3 * 9 * nn + 9 * dd] +=
840  g_w[gg] * N_face[nb_dofs_face * gg + nn] * s2_imag[dd] *
841  traction[1];
842  }
843  }
844  }
845  }
846  }
847  MoFEMFunctionReturnHot(0);
848 }
849 MoFEMErrorCode KExt_hh_hierarchical_edge(
850  double eps, int order, int *order_edge, double *N, double *N_face,
851  double *N_edge[], double *diffN, double *diffN_face, double *diffN_edge[],
852  double *t, double *t_edge[], double *t_face, double *dofs_x,
853  double *dofs_x_edge[], double *dofs_x_face, double *KExt_hedge[3],
854  double *KExt_edgeedge[3][3], double *KExt_faceedge[3], int g_dim,
855  const double *g_w) {
856  MoFEMFunctionBeginHot;
857  int gg, dd, ii, nn, ee, EE;
858  int nb_dofs_face = NBFACETRI_H1(order);
859  for (EE = 0; EE < 3; EE++) {
860  int nb_dofs_edge_EE = NBEDGE_H1(order_edge[EE]);
861  bzero(KExt_hedge[EE], 9 * 3 * nb_dofs_edge_EE * sizeof(double));
862  if (KExt_edgeedge != NULL) {
863  for (ee = 0; ee < 3; ee++) {
864  int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
865  bzero(KExt_edgeedge[EE][ee],
866  3 * nb_dofs_edge_EE * 3 * nb_dofs_edge * sizeof(double));
867  }
868  }
869  if (KExt_faceedge != NULL) {
870  bzero(KExt_faceedge[EE],
871  3 * nb_dofs_edge_EE * 3 * nb_dofs_face * sizeof(double));
872  }
873  }
874  for (gg = 0; gg < g_dim; gg++) {
875  double traction[3] = {0, 0, 0};
876  CHKERR Traction_hierarchical(order, order_edge, N, N_face, N_edge, t,
877  t_edge, t_face, traction, gg);
878  for (EE = 0; EE < 3; EE++) {
879  int nb_dofs_edge_EE = NBEDGE_H1(order_edge[EE]);
880  double *idofs_x_edge[3] = {NULL, NULL, NULL};
881  double idofs_x_edge_EE[3 * nb_dofs_edge_EE];
882  idofs_x_edge[EE] = idofs_x_edge_EE;
883  for (ii = 0; ii < 3 * nb_dofs_edge_EE; ii++) {
884  bzero(idofs_x_edge_EE, 3 * nb_dofs_edge_EE * sizeof(double));
885  idofs_x_edge_EE[ii] = eps;
886  __CLPK_doublecomplex xnormal[3], xs1[3], xs2[3];
887  CHKERR Normal_hierarchical(order, order_edge, // FIXME
888  order, order_edge, diffN, diffN_face,
889  diffN_edge, dofs_x, dofs_x_edge, dofs_x_face,
890  NULL, idofs_x_edge, NULL, xnormal, xs1, xs2,
891  gg);
892  CHKERR Base_scale(xnormal, xs1, xs2);
893  double normal_imag[3], s1_imag[3], s2_imag[3];
894  for (dd = 0; dd < 3; dd++) {
895  normal_imag[dd] = 0.5 * xnormal[dd].i / eps;
896  s1_imag[dd] = 0.5 * xs1[dd].i / eps;
897  s2_imag[dd] = 0.5 * xs2[dd].i / eps;
898  }
899  for (nn = 0; nn < 3; nn++) {
900  for (dd = 0; dd < 3; dd++) {
901  KExt_hedge[EE][ii + 3 * nb_dofs_edge_EE * 3 * nn +
902  3 * nb_dofs_edge_EE * dd] +=
903  g_w[gg] * N[3 * gg + nn] * normal_imag[dd] * traction[2];
904  KExt_hedge[EE][ii + 3 * nb_dofs_edge_EE * 3 * nn +
905  3 * nb_dofs_edge_EE * dd] +=
906  g_w[gg] * N[3 * gg + nn] * s1_imag[dd] * traction[0];
907  KExt_hedge[EE][ii + 3 * nb_dofs_edge_EE * 3 * nn +
908  3 * nb_dofs_edge_EE * dd] +=
909  g_w[gg] * N[3 * gg + nn] * s2_imag[dd] * traction[1];
910  }
911  }
912  if (KExt_edgeedge != NULL) {
913  for (ee = 0; ee < 3; ee++) {
914  int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
915  for (nn = 0; nn < nb_dofs_edge; nn++) {
916  for (dd = 0; dd < 3; dd++) {
917  KExt_edgeedge[EE][ee][ii + 3 * nb_dofs_edge_EE * 3 * nn +
918  3 * nb_dofs_edge_EE * dd] +=
919  g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] *
920  normal_imag[dd] * traction[2];
921  KExt_edgeedge[EE][ee][ii + 3 * nb_dofs_edge_EE * 3 * nn +
922  3 * nb_dofs_edge_EE * dd] +=
923  g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] * s1_imag[dd] *
924  traction[0];
925  KExt_edgeedge[EE][ee][ii + 3 * nb_dofs_edge_EE * 3 * nn +
926  3 * nb_dofs_edge_EE * dd] +=
927  g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] * s2_imag[dd] *
928  traction[1];
929  }
930  }
931  }
932  }
933  if (KExt_faceedge != NULL) {
934  for (nn = 0; nn < nb_dofs_face; nn++) {
935  for (dd = 0; dd < 3; dd++) {
936  KExt_faceedge[EE][ii + 3 * nb_dofs_edge_EE * 3 * nn +
937  3 * nb_dofs_edge_EE * dd] +=
938  g_w[gg] * N_face[nb_dofs_face * gg + nn] * normal_imag[dd] *
939  traction[2];
940  KExt_faceedge[EE][ii + 3 * nb_dofs_edge_EE * 3 * nn +
941  3 * nb_dofs_edge_EE * dd] +=
942  g_w[gg] * N_face[nb_dofs_face * gg + nn] * s1_imag[dd] *
943  traction[0];
944  KExt_faceedge[EE][ii + 3 * nb_dofs_edge_EE * 3 * nn +
945  3 * nb_dofs_edge_EE * dd] +=
946  g_w[gg] * N_face[nb_dofs_face * gg + nn] * s2_imag[dd] *
947  traction[1];
948  }
949  }
950  }
951  }
952  }
953  }
954  MoFEMFunctionReturnHot(0);
955 }
956 MoFEMErrorCode
957 KExt_hh_hierarchical_face(double eps, int order, int *order_edge, double *N,
958  double *N_face, double *N_edge[], double *diffN,
959  double *diffN_face, double *diffN_edge[], double *t,
960  double *t_edge[], double *t_face, double *dofs_x,
961  double *dofs_x_edge[], double *dofs_x_face,
962  double *KExt_hface, double *KExt_edgeface[3],
963  double *KExt_faceface, int g_dim, const double *g_w) {
964  MoFEMFunctionBeginHot;
965  int gg, dd, ii, nn, ee;
966  int nb_dofs_face = NBFACETRI_H1(order);
967  bzero(KExt_hface, 9 * 3 * nb_dofs_face * sizeof(double));
968  if (KExt_edgeface != NULL) {
969  for (ee = 0; ee < 3; ee++) {
970  int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
971  bzero(KExt_edgeface[ee],
972  3 * nb_dofs_face * 3 * nb_dofs_edge * sizeof(double));
973  }
974  }
975  if (KExt_faceface != NULL) {
976  bzero(KExt_faceface, 3 * nb_dofs_face * 3 * nb_dofs_face * sizeof(double));
977  }
978  for (gg = 0; gg < g_dim; gg++) {
979  double traction[3] = {0, 0, 0};
980  CHKERR Traction_hierarchical(order, order_edge, N, N_face, N_edge, t,
981  t_edge, t_face, traction, gg);
982  double idofs_x_face[3 * nb_dofs_face];
983  for (ii = 0; ii < 3 * nb_dofs_face; ii++) {
984  bzero(idofs_x_face, 3 * nb_dofs_face * sizeof(double));
985  idofs_x_face[ii] = eps;
986  __CLPK_doublecomplex xnormal[3], xs1[3], xs2[3];
987  CHKERR Normal_hierarchical(
988  order, order_edge, // FIXME
989  order, order_edge, diffN, diffN_face, diffN_edge, dofs_x, dofs_x_edge,
990  dofs_x_face, NULL, NULL, idofs_x_face, xnormal, xs1, xs2, gg);
991  CHKERR Base_scale(xnormal, xs1, xs2);
992  double normal_imag[3], s1_imag[3], s2_imag[3];
993  for (dd = 0; dd < 3; dd++) {
994  normal_imag[dd] = 0.5 * xnormal[dd].i / eps;
995  s1_imag[dd] = 0.5 * xs1[dd].i / eps;
996  s2_imag[dd] = 0.5 * xs2[dd].i / eps;
997  }
998  for (nn = 0; nn < 3; nn++) {
999  for (dd = 0; dd < 3; dd++) {
1000  KExt_hface[ii + 3 * nb_dofs_face * 3 * nn + 3 * nb_dofs_face * dd] +=
1001  g_w[gg] * N[3 * gg + nn] * normal_imag[dd] * traction[2];
1002  KExt_hface[ii + 3 * nb_dofs_face * 3 * nn + 3 * nb_dofs_face * dd] +=
1003  g_w[gg] * N[3 * gg + nn] * s1_imag[dd] * traction[0];
1004  KExt_hface[ii + 3 * nb_dofs_face * 3 * nn + 3 * nb_dofs_face * dd] +=
1005  g_w[gg] * N[3 * gg + nn] * s2_imag[dd] * traction[1];
1006  }
1007  }
1008  if (KExt_edgeface != NULL) {
1009  for (ee = 0; ee < 3; ee++) {
1010  int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
1011  for (nn = 0; nn < nb_dofs_edge; nn++) {
1012  for (dd = 0; dd < 3; dd++) {
1013  KExt_edgeface[ee][ii + 3 * nb_dofs_face * 3 * nn +
1014  3 * nb_dofs_face * dd] +=
1015  g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] *
1016  normal_imag[dd] * traction[2];
1017  KExt_edgeface[ee][ii + 3 * nb_dofs_face * 3 * nn +
1018  3 * nb_dofs_face * dd] +=
1019  g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] * s1_imag[dd] *
1020  traction[0];
1021  KExt_edgeface[ee][ii + 3 * nb_dofs_face * 3 * nn +
1022  3 * nb_dofs_face * dd] +=
1023  g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] * s2_imag[dd] *
1024  traction[1];
1025  }
1026  }
1027  }
1028  }
1029  if (KExt_faceface != NULL) {
1030  for (nn = 0; nn < nb_dofs_face; nn++) {
1031  for (dd = 0; dd < 3; dd++) {
1032  KExt_faceface[ii + 3 * nb_dofs_face * 3 * nn +
1033  3 * nb_dofs_face * dd] +=
1034  g_w[gg] * N_face[nb_dofs_face * gg + nn] * normal_imag[dd] *
1035  traction[2];
1036  KExt_faceface[ii + 3 * nb_dofs_face * 3 * nn +
1037  3 * nb_dofs_face * dd] +=
1038  g_w[gg] * N_face[nb_dofs_face * gg + nn] * s1_imag[dd] *
1039  traction[0];
1040  KExt_faceface[ii + 3 * nb_dofs_face * 3 * nn +
1041  3 * nb_dofs_face * dd] +=
1042  g_w[gg] * N_face[nb_dofs_face * gg + nn] * s2_imag[dd] *
1043  traction[1];
1044  }
1045  }
1046  }
1047  }
1048  }
1049  MoFEMFunctionReturnHot(0);
1050 }
1051 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::lHs
virtual MoFEMErrorCode lHs()
Definition: SurfacePressureComplexForLazy.cpp:264
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::calcTraction
virtual MoFEMErrorCode calcTraction()
Definition: SurfacePressureComplexForLazy.cpp:396
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
SIDESET
@ SIDESET
Definition: definitions.h:160
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: SurfacePressureComplexForLazy.cpp:443
NeumannForcesSurfaceComplexForLazy::spatialField
const std::string spatialField
Definition: SurfacePressureComplexForLazy.hpp:175
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1644
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::MyTriangleSpatialFE
MyTriangleSpatialFE(MoFEM::Interface &_mField, Mat _Aij, Vec &_F, double *scale_lhs, double *scale_rhs, std::string spatial_field_name="SPATIAL_POSITION", std::string mat_field_name="MESH_NODE_POSITIONS")
Definition: SurfacePressureComplexForLazy.cpp:161
NBEDGE_H1
#define NBEDGE_H1(P)
Number of base function on edge for H1 space.
Definition: h1_hdiv_hcurl_l2.h:55
EntityHandle
MoFEM::EntitiesFieldData::EntData::getOrder
ApproximationOrder getOrder() const
get approximation order
Definition: EntitiesFieldData.hpp:1210
NeumannForcesSurfaceComplexForLazy::sCale
double * sCale
Definition: SurfacePressureComplexForLazy.hpp:156
MoFEM::FaceElementForcesAndSourcesCore::meshPositionsFieldName
std::string meshPositionsFieldName
Definition: FaceElementForcesAndSourcesCore.hpp:29
NeumannForcesSurfaceComplexForLazy::AuxMethodMaterial
Definition: SurfacePressureComplexForLazy.hpp:29
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::sCaleLhs
double * sCaleLhs
Definition: SurfacePressureComplexForLazy.hpp:41
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::rHs
virtual MoFEMErrorCode rHs()
Definition: SurfacePressureComplexForLazy.cpp:185
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
_F
#define _F(n)
Definition: quad.c:25
lapack_dgesv
static __CLPK_integer lapack_dgesv(__CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb)
Definition: lapack_wrap.h:176
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
Base_scale
PetscErrorCode Base_scale(__CLPK_doublecomplex *xnormal, __CLPK_doublecomplex *xs1, __CLPK_doublecomplex *xs2)
Definition: fem_tools.c:741
MoFEM::CubitMeshSets
this struct keeps basic methods for moab meshset about material and boundary conditions
Definition: BCMultiIndices.hpp:19
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::SnesMethod::snes_B
Mat & snes_B
preconditioner of jacobian matrix
Definition: LoopMethods.hpp:124
__CLPK_doublecomplex::r
__CLPK_doublereal r
Definition: lapack_wrap.h:35
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::addForce
MoFEMErrorCode addForce(int ms_id)
Definition: SurfacePressureComplexForLazy.cpp:516
MoFEM::EntitiesFieldData::EntData::getFieldData
const VectorDouble & getFieldData() const
get dofs values
Definition: EntitiesFieldData.hpp:1254
NeumannForcesSurfaceComplexForLazy::materialField
const std::string materialField
Definition: SurfacePressureComplexForLazy.hpp:176
MoFEM.hpp
MOAB_THROW
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:554
MoFEM::EntitiesFieldData::EntData::getDiffN
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
Definition: EntitiesFieldData.hpp:1329
tricircumcenter3d_tp
void tricircumcenter3d_tp(double a[3], double b[3], double c[3], double circumcenter[3], double *xi, double *eta)
NeumannForcesSurfaceComplexForLazy::AuxMethodMaterial::AuxMethodMaterial
AuxMethodMaterial(const string &field_name, MyTriangleSpatialFE *my_ptr, const char type)
Definition: SurfacePressureComplexForLazy.cpp:61
scale
double scale
Definition: plastic.cpp:119
tetcircumcenter_tp
void tetcircumcenter_tp(double a[3], double b[3], double c[3], double d[3], double circumcenter[3], double *xi, double *eta, double *zeta)
NeumannForcesSurfaceComplexForLazy::NeumannForcesSurfaceComplexForLazy
NeumannForcesSurfaceComplexForLazy(MoFEM::Interface &m_field, Mat _Aij, Vec _F, double *scale_lhs, double *scale_rhs, std::string spatial_field_name="SPATIAL_POSITION", std::string material_field_name="MESH_NODE_POSITIONS")
Definition: SurfacePressureComplexForLazy.cpp:576
zeta
double zeta
Viscous hardening.
Definition: plastic.cpp:126
MoFEM::CubitMeshSets::getBcDataStructure
MoFEMErrorCode getBcDataStructure(CUBIT_BC_DATA_TYPE &data) const
Definition: BCMultiIndices.hpp:296
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
NeumannForcesSurfaceComplexForLazy::feSpatial
MyTriangleSpatialFE feSpatial
Definition: SurfacePressureComplexForLazy.hpp:152
sdf.r
int r
Definition: sdf.py:8
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
KExt_hh_hierarchical
MoFEMErrorCode KExt_hh_hierarchical(double eps, int order, int *order_edge, double *N, double *N_face, double *N_edge[], double *diffN, double *diffN_face, double *diffN_edge[], double *t, double *t_edge[], double *t_face, double *dofs_x, double *dofs_x_edge[], double *dofs_x_face, double *KExt_hh, double *KExt_egdeh[3], double *KExt_faceh, int g_dim, const double *g_w)
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
NODESET
@ NODESET
Definition: definitions.h:159
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::F
Vec F
Definition: SurfacePressureComplexForLazy.hpp:51
Traction_hierarchical
MoFEMErrorCode Traction_hierarchical(int order, int *order_edge, double *N, double *N_face, double *N_edge[], double *t, double *t_edge[], double *t_face, double *traction, int gg)
Definition: SurfacePressureComplexForLazy.cpp:587
SurfacePressureComplexForLazy.hpp
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
NeumannForcesSurfaceComplexForLazy::AuxMethodSpatial
Definition: SurfacePressureComplexForLazy.hpp:19
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::sCaleRhs
double * sCaleRhs
Definition: SurfacePressureComplexForLazy.hpp:42
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
KExt_hh_hierarchical_edge
MoFEMErrorCode KExt_hh_hierarchical_edge(double eps, int order, int *order_edge, double *N, double *N_face, double *N_edge[], double *diffN, double *diffN_face, double *diffN_edge[], double *t, double *t_edge[], double *t_face, double *dofs_x, double *dofs_x_edge[], double *dofs_x_face, double *Khext_edge[3], double *KExt_edgeegde[3][3], double *KExt_faceedge[3], int g_dim, const double *g_w)
Definition: SurfacePressureComplexForLazy.cpp:849
convert.type
type
Definition: convert.py:64
Fext_h_hierarchical
MoFEMErrorCode Fext_h_hierarchical(int order, int *order_edge, double *N, double *N_face, double *N_edge[], double *diffN, double *diffN_face, double *diffN_edge[], double *t, double *t_edge[], double *t_face, double *dofs_x, double *dofs_x_edge[], double *dofs_x_face, double *idofs_x, double *idofs_x_edge[], double *idofs_x_face, double *Fext, double *Fext_egde[], double *Fext_face, double *iFext, double *iFext_egde[], double *iFext_face, int g_dim, const double *g_w)
Definition: SurfacePressureComplexForLazy.cpp:620
NeumannForcesSurfaceComplexForLazy::AuxMethodMaterial::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: SurfacePressureComplexForLazy.cpp:117
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1214
eta
double eta
Definition: free_surface.cpp:170
NeumannForcesSurfaceComplexForLazy::AuxMethodSpatial::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: SurfacePressureComplexForLazy.cpp:66
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MoFEM::CubitMeshSets::meshset
EntityHandle meshset
Definition: BCMultiIndices.hpp:21
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE
Definition: SurfacePressureComplexForLazy.hpp:39
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::operator()
MoFEMErrorCode operator()()
function is run for every finite element
Definition: SurfacePressureComplexForLazy.cpp:461
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::methodsOp
boost::ptr_vector< MethodForForceScaling > methodsOp
Definition: SurfacePressureComplexForLazy.hpp:148
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
KExt_hh_hierarchical_face
MoFEMErrorCode KExt_hh_hierarchical_face(double eps, int order, int *order_edge, double *N, double *N_face, double *N_edge[], double *diffN, double *diffN_face, double *diffN_edge[], double *t, double *t_edge[], double *t_face, double *dofs_x, double *dofs_x_edge[], double *dofs_x_face, double *KExt_hface, double *KExt_egdeface[3], double *KExt_faceface, int g_dim, const double *g_w)
Definition: SurfacePressureComplexForLazy.cpp:957
N
const int N
Definition: speed_test.cpp:3
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::Aij
Mat Aij
Definition: SurfacePressureComplexForLazy.hpp:50
MoFEM::SnesMethod::snes_f
Vec & snes_f
residual
Definition: LoopMethods.hpp:122
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
NBFACETRI_H1
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
Definition: h1_hdiv_hcurl_l2.h:60
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: EntitiesFieldData.hpp:1318
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:83
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
ShapeFaceBaseMBTRI
PetscErrorCode ShapeFaceBaseMBTRI(double *diffN, const double *coords, double *normal, double *s1, double *s2)
Definition: fem_tools.c:204
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::addPressure
MoFEMErrorCode addPressure(int ms_id)
Definition: SurfacePressureComplexForLazy.cpp:531
Normal_hierarchical
PetscErrorCode Normal_hierarchical(int order_approx, int *order_edge_approx, int order, int *order_edge, double *diffN, double *diffN_face, double *diffN_edge[], double *dofs, double *dofs_edge[], double *dofs_face, double *idofs, double *idofs_edge[], double *idofs_face, __CLPK_doublecomplex *xnormal, __CLPK_doublecomplex *s1, __CLPK_doublecomplex *s2, int gg)
Complex normal.
Definition: fem_tools.c:569
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
NeumannForcesSurfaceComplexForLazy::thScale
Tag thScale
Definition: SurfacePressureComplexForLazy.hpp:154
NeumannForcesSurfaceComplexForLazy::MyTriangleSpatialFE::reBaseToFaceLoocalCoordSystem
MoFEMErrorCode reBaseToFaceLoocalCoordSystem(MatrixDouble &t_glob_nodal)
Definition: SurfacePressureComplexForLazy.cpp:377
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
__CLPK_integer
long int __CLPK_integer
Definition: lapack_wrap.h:23
MethodForForceScaling::applyScale
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
Definition: MethodForForceScaling.hpp:21
MoFEM::MeshsetsManager::getCubitMeshsetPtr
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
Definition: MeshsetsManager.cpp:578
__CLPK_doublecomplex::i
__CLPK_doublereal i
Definition: lapack_wrap.h:35
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
__CLPK_doublecomplex
Definition: lapack_wrap.h:34
NeumannForcesSurfaceComplexForLazy::AuxMethodSpatial::AuxMethodSpatial
AuxMethodSpatial(const string &field_name, MyTriangleSpatialFE *my_ptr, const char type)
Definition: SurfacePressureComplexForLazy.cpp:56
SurfacePressure.hpp
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
NeumannForcesSurfaceComplexForLazy::mField
MoFEM::Interface & mField
Definition: SurfacePressureComplexForLazy.hpp:151
MoFEM::ForcesAndSourcesCore::mField
Interface & mField
Definition: ForcesAndSourcesCore.hpp:24
F
@ F
Definition: free_surface.cpp:394
Generated by Doxygen 1.8.17 and hosted at University of Glasgow