v0.14.0
MortarContactProblem.cpp
Go to the documentation of this file.
1 /* This file is part of MoFEM.
2  * MoFEM is free software: you can redistribute it and/or modify it under
3  * the terms of the GNU Lesser General Public License as published by the
4  * Free Software Foundation, either version 3 of the License, or (at your
5  * option) any later version.
6  *
7  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
8  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
10  * License for more details.
11  *
12  * You should have received a copy of the GNU Lesser General Public
13  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
14 
15 MoFEMErrorCode
16 MortarContactProblem::MortarContactElement::setGaussPts(int order) {
17  MoFEMFunctionBegin;
18 
19  FTensor::Index<'i', 3> i;
20  FTensor::Index<'j', 3> j;
21  auto get_tensor_vec = [](VectorDouble &v) {
22  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v(0), &v(1),
23  &v(2));
24  };
25 
26  order *= 2; // multiply by 2 due to the integrand of NTN (twice the
27  // approximation)
28 
29  // Defined integration points for only 1 integrated tris
30  int nb_gauss_pts_1tri;
31  MatrixDouble gaussPts_1tri;
32  nb_gauss_pts_1tri = QUAD_2D_TABLE[order]->npoints;
33 
34  gaussPts_1tri.resize(3, nb_gauss_pts_1tri, false);
35  cblas_dcopy(nb_gauss_pts_1tri, &QUAD_2D_TABLE[order]->points[1], 3,
36  &gaussPts_1tri(0, 0), 1);
37  cblas_dcopy(nb_gauss_pts_1tri, &QUAD_2D_TABLE[order]->points[2], 3,
38  &gaussPts_1tri(1, 0), 1);
39  cblas_dcopy(nb_gauss_pts_1tri, QUAD_2D_TABLE[order]->weights, 1,
40  &gaussPts_1tri(2, 0), 1);
41  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts_1tri, 3,
42  false);
43  double *shape_ptr =
44  &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
45  cblas_dcopy(3 * nb_gauss_pts_1tri, QUAD_2D_TABLE[order]->points, 1, shape_ptr,
46  1);
47 
48  // get the entity (prism in this case)
49  EntityHandle prism_entity =
50  numeredEntFiniteElementPtr->getEnt(); // we will get prism element here
51  typedef ContactSearchKdTree::ContactCommonData_multiIndex::index<
52  ContactSearchKdTree::Prism_tag>::type::iterator ItMultIndexPrism;
53  ItMultIndexPrism it_mult_index_prism =
54  contactCommondataMultiIndex->get<ContactSearchKdTree::Prism_tag>().find(
55  prism_entity);
56 
57  // get integration tris of each prism
58  Range range_poly_tris;
59  range_poly_tris.clear();
60  range_poly_tris = it_mult_index_prism->get()->commonIntegratedTriangle;
61 
62  // Number of integration points = number of Gauss points in each common
63  // tris * no of common tris * 2 (2 here is due to different integratio
64  // rule
65  // of master and slave tris)
66  int nb_gauss_pts = nb_gauss_pts_1tri * range_poly_tris.size();
67  gaussPtsMaster.resize(3, nb_gauss_pts, false);
68  gaussPtsSlave.resize(3, nb_gauss_pts, false);
69 
70  const EntityHandle *conn_slave = NULL;
71  int num_nodes_prism = 0;
72  rval = mField.get_moab().get_connectivity(prism_entity, conn_slave,
73  num_nodes_prism);
74  CHKERRQ_MOAB(rval);
75  // for prism element bottom 3 nodes belong to slave tri and top 3 nodes
76  // belong to master tri
77  VectorDouble coords_prism; // 6*3=18 coordinates for 6 nodes in 3D
78  coords_prism.resize(18, false);
79  rval = mField.get_moab().get_coords(conn_slave, num_nodes_prism,
80  &*coords_prism.data().begin());
81  CHKERRQ_MOAB(rval);
82 
83  VectorDouble v_elem_coords_master, v_elem_coords_slave;
84  v_elem_coords_master.resize(9, false);
85  v_elem_coords_master = subrange(coords_prism, 0, 9);
86  v_elem_coords_slave.resize(9, false);
87  v_elem_coords_slave = subrange(coords_prism, 9, 18);
88  auto t_elem_coords_master = get_tensor_vec(v_elem_coords_master);
89  auto t_elem_coords_slave = get_tensor_vec(v_elem_coords_slave);
90  // Here we need to calculate the local ara of integration triangles in
91  // both master and slave surfaces (we need this for weight calculation of
92  // the Gauss points)
93 
94  // To do this, first convert the master and slave triangles to 2D case
95  // (or z=0) as both master and slave triangles are oriented in 3D space
96 
97  // tansfermation(rotation) matrix (use the same matrix as defined in
98  // ContactSearchKdTree)
99  MatrixDouble m_rot;
100  m_rot.resize(3, 3, false);
101 
102  FTensor::Tensor2<double *, 3, 3> t_m_rot{
103  &m_rot(0, 0), &m_rot(0, 1), &m_rot(0, 2), &m_rot(1, 0), &m_rot(1, 1),
104  &m_rot(1, 2), &m_rot(2, 0), &m_rot(2, 1), &m_rot(2, 2)};
105 
106  ContactSearchKdTree contact_problem_mulit_index(mField);
107 
108  CHKERR contact_problem_mulit_index.rotationMatrix(m_rot, v_elem_coords_slave);
109 
110  VectorDouble v_elem_coords_master_new, v_elem_coords_slave_new;
111  v_elem_coords_master_new.resize(9, false);
112  v_elem_coords_slave_new.resize(9, false);
113  auto t_elem_coords_master_new = get_tensor_vec(v_elem_coords_master_new);
114  auto t_elem_coords_slave_new = get_tensor_vec(v_elem_coords_slave_new);
115 
116  // master and slave tri elemnets (coord in 2D)
117  VectorDouble v_elem_coords_master_2d, v_elem_coords_slave_2d;
118  v_elem_coords_master_2d.resize(6, false);
119  v_elem_coords_slave_2d.resize(6, false);
120  FTensor::Tensor1<double *, 2> t_elem_coords_master_2d{
121  &v_elem_coords_master_2d(0), &v_elem_coords_master_2d(1), 2};
122 
123  FTensor::Tensor1<double *, 2> t_elem_coords_slave_2d{
124  &v_elem_coords_slave_2d(0), &v_elem_coords_slave_2d(1), 2};
125  // transfer the slave and master tri coordinates (n=3) one by one to z
126  // plane
127  for (int ii = 0; ii != 3; ++ii) {
128  t_elem_coords_slave_new(i) = t_m_rot(i, j) * t_elem_coords_slave(j);
129  t_elem_coords_master_new(i) = t_m_rot(i, j) * t_elem_coords_master(j);
130 
131  t_elem_coords_master_2d(0) = t_elem_coords_master_new(0);
132  t_elem_coords_master_2d(1) = t_elem_coords_master_new(1);
133 
134  t_elem_coords_slave_2d(0) = t_elem_coords_slave_new(0);
135  t_elem_coords_slave_2d(1) = t_elem_coords_slave_new(1);
136 
137  ++t_elem_coords_slave_new;
138  ++t_elem_coords_slave;
139  ++t_elem_coords_master_new;
140  ++t_elem_coords_master;
141 
142  ++t_elem_coords_master_2d;
143  ++t_elem_coords_slave_2d;
144  }
145 
146  // get nodal coordinates
147  VectorDouble v_coords_integration_tri;
148  v_coords_integration_tri.resize(9, false);
149 
150  // transfer coord to z plane
151  VectorDouble v_coords_integration_tri_new;
152  v_coords_integration_tri_new.resize(9, false);
153 
154  VectorDouble loc_coords_master; // starting point
155  loc_coords_master.resize(2, false);
156  double n_input[3] = {1, 0, 0}; // shape function at starting point
157 
158  // function to calculate the local coordinate of element based on its
159  // global coordinates master is the bottom surface
160  VectorDouble loc_coords_slave;
161  loc_coords_slave.resize(2, false);
162 
163  // global coordinate of each Gauss point in 2D (z plane)
164  // x and y global coordinates
165  VectorDouble v_glob_coords;
166  v_glob_coords.resize(2, false);
167 
168  auto t_gaussPtsMaster = getFTensor1FromMat<3>(gaussPtsMaster);
169  auto t_gaussPtsSlave = getFTensor1FromMat<3>(gaussPtsSlave);
170  // for each prism loop over all the integration tris
171  for (Range::iterator it_tri = range_poly_tris.begin();
172  it_tri != range_poly_tris.end(); ++it_tri) {
173  const EntityHandle *conn_face;
174  int num_nodes_tri;
175  // get nodes attached to the tri
176  rval = mField.get_moab().get_connectivity(*it_tri, conn_face, num_nodes_tri,
177  true);
178  CHKERRQ_MOAB(rval);
179  // v_coords_integration_tri.clear(); //[x1 y1 z1 x2 y2 z2 .......]
180  rval = mField.get_moab().get_coords(
181  conn_face, num_nodes_tri, &*v_coords_integration_tri.data().begin());
182  CHKERRQ_MOAB(rval);
183 
184  auto t_coords_integration_tri = get_tensor_vec(v_coords_integration_tri);
185  auto t_coords_integration_tri_new =
186  get_tensor_vec(v_coords_integration_tri_new);
187 
188  for (int ii = 0; ii != 3; ++ii) {
189  t_coords_integration_tri_new(i) =
190  t_m_rot(i, j) * t_coords_integration_tri(j);
191  t_coords_integration_tri_new(2) = 0;
192  ++t_coords_integration_tri;
193  ++t_coords_integration_tri_new;
194  }
195 
196  // shape function derivative for tri elements (these are constant)
197  // diff_n_tri = [dN1/dxi, dN1/deta, dN2/dxi, dN2/deta, dN3/dxi,
198  // dN3/deta] = [-1 -1 1 0 0 1]
199  double diff_n_tri[6];
200  CHKERR ShapeDiffMBTRI(diff_n_tri);
201 
202  // calculate local coordinates of each integration tri
203  // double n_input[3] = {1, 0, 0}; // shape function at starting point
204  // function to calculate the local coordinate of element based on its
205  // global coordinates element local coordinates of nodes of integration
206  // tri
207  // double coords_integration_tri_loc_master[9],
208  // coords_integration_tri_loc_slave[9];
209 
210  VectorDouble coords_integration_tri_loc_master,
211  coords_integration_tri_loc_slave;
212  coords_integration_tri_loc_master.resize(6, false);
213  coords_integration_tri_loc_slave.resize(6, false);
214  FTensor::Tensor1<double *, 2> t_coords_integration_tri_loc_master{
215  &coords_integration_tri_loc_master(0),
216  &coords_integration_tri_loc_master(1), 2};
217  FTensor::Tensor1<double *, 2> t_coords_integration_tri_loc_slave{
218  &coords_integration_tri_loc_slave(0),
219  &coords_integration_tri_loc_slave(1), 2};
220 
221  auto t_coords_integration_tri_new_2 =
222  get_tensor_vec(v_coords_integration_tri_new);
223 
224  for (int ii = 0; ii != 3; ++ii) {
225 
226  double glob_coords_tri[2];
227  glob_coords_tri[0] = t_coords_integration_tri_new_2(0);
228  glob_coords_tri[1] = t_coords_integration_tri_new_2(1);
229 
230  // local coordinates of integration tri in master element
231  double loc_coords_tri_master[2] = {0, 0}; // starting point
232  CHKERR ShapeMBTRI_inverse(n_input, diff_n_tri,
233  &v_elem_coords_master_2d(0), glob_coords_tri,
234  loc_coords_tri_master);
235 
236  t_coords_integration_tri_loc_master(0) = loc_coords_tri_master[0];
237  t_coords_integration_tri_loc_master(1) = loc_coords_tri_master[1];
238 
239  // local coordinates of integration tri in slave element
240  double loc_coords_tri_slave[2] = {0, 0}; // starting point
241  CHKERR ShapeMBTRI_inverse(n_input, diff_n_tri, &v_elem_coords_slave_2d(0),
242  glob_coords_tri, loc_coords_tri_slave);
243 
244  t_coords_integration_tri_loc_slave(0) = loc_coords_tri_slave[0];
245  t_coords_integration_tri_loc_slave(1) = loc_coords_tri_slave[1];
246 
247  ++t_coords_integration_tri_loc_master;
248  ++t_coords_integration_tri_loc_slave;
249  ++t_coords_integration_tri_new_2;
250  }
251 
252  // local (not global) area or jacobian of integration tri in both
253  // master and slave triangles (can be +ve or -ve based on surface
254  // orientation)
255  double area_integration_tri_master_loc, area_integration_tri_slave_loc;
256  area_integration_tri_master_loc =
257  std::abs(area2D(&coords_integration_tri_loc_master(0),
258  &coords_integration_tri_loc_master(2),
259  &coords_integration_tri_loc_master(4)));
260 
261  area_integration_tri_slave_loc =
262  std::abs(area2D(&coords_integration_tri_loc_slave(0),
263  &coords_integration_tri_loc_slave(2),
264  &coords_integration_tri_loc_slave(4)));
265 
266  auto t_gaussPts_1tri = getFTensor1FromMat<3>(gaussPts_1tri);
267 
268  // for each integration tri loop over all its Gauss points
269  // calculate global coordinates of each integration point and then
270  // calculate the local coordinates of this integration point in each
271  // master and slave surface
272  for (int gg = 0; gg != nb_gauss_pts_1tri; ++gg) {
273 
274  t_gaussPtsMaster(2) =
275  t_gaussPts_1tri(2) * area_integration_tri_master_loc *
276  2; // 2 here is due to as for ref tri A=1/2 and w=1 or w=2*A
277  t_gaussPtsSlave(2) =
278  t_gaussPts_1tri(2) * area_integration_tri_slave_loc * 2;
279 
280  // shape function for each Gauss point
281  MatrixDouble N_tri;
282  N_tri.resize(1, 3, false);
283 
284  CHKERR ShapeMBTRI(&N_tri(0, 0), &t_gaussPts_1tri(0), &t_gaussPts_1tri(1),
285  1);
286 
287  FTensor::Tensor1<double *, 3> t_N_tri{&N_tri(0, 0), &N_tri(0, 1),
288  &N_tri(0, 2)};
289 
290  FTensor::Tensor1<double *, 3> t_coords_integration_tri_new{
291  &v_coords_integration_tri_new(0), &v_coords_integration_tri_new(3),
292  &v_coords_integration_tri_new(6), 1};
293 
294  v_glob_coords(0) = t_N_tri(i) * t_coords_integration_tri_new(i);
295  ++t_coords_integration_tri_new;
296  v_glob_coords(1) = t_N_tri(i) * t_coords_integration_tri_new(i);
297 
298  // calculate local coordinates of each Gauss point in both slave
299  // and master tris
300  CHKERR ShapeMBTRI_inverse(
301  n_input, diff_n_tri, &v_elem_coords_master_2d(0),
302  &*v_glob_coords.data().begin(), &loc_coords_master(0));
303 
304  t_gaussPtsMaster(0) = loc_coords_master(0);
305  t_gaussPtsMaster(1) = loc_coords_master(1);
306 
307  // slave is the top surface
308  CHKERR ShapeMBTRI_inverse(n_input, diff_n_tri, &v_elem_coords_slave_2d(0),
309  &*v_glob_coords.data().begin(),
310  &loc_coords_slave(0));
311 
312  t_gaussPtsSlave(0) = loc_coords_slave[0];
313  t_gaussPtsSlave(1) = loc_coords_slave[1];
314 
315  ++t_gaussPts_1tri;
316  ++t_gaussPtsSlave;
317  ++t_gaussPtsMaster;
318  }
319  }
320  MoFEMFunctionReturn(0);
321 }
322 
323 MoFEMErrorCode
324 MortarContactProblem::MortarConvectMasterContactElement::setGaussPts(
325  int order) {
326  MoFEMFunctionBegin;
327  CHKERR MortarContactElement::setGaussPts(order);
328  CHKERR convectPtr->convectSlaveIntegrationPts<true>();
329  MoFEMFunctionReturn(0);
330 }
331 
332 MoFEMErrorCode
333 MortarContactProblem::MortarConvectSlaveContactElement::setGaussPts(int order) {
334  MoFEMFunctionBegin;
335  CHKERR MortarContactElement::setGaussPts(order);
336  CHKERR convectPtr->convectSlaveIntegrationPts<false>();
337  MoFEMFunctionReturn(0);
338 }
339 
340 MoFEMErrorCode MortarContactProblem::setContactOperatorsRhs(
341  boost::shared_ptr<MortarContactElement> fe_rhs_mortar_contact,
342  boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
343  string field_name, string lagrange_field_name, bool is_alm,
344  bool is_eigen_pos_field, string eigen_pos_field_name,
345  bool use_reference_coordinates) {
346  MoFEMFunctionBegin;
347 
348  fe_rhs_mortar_contact->getOpPtrVector().push_back(new OpGetNormalSlaveALE(
349  "MESH_NODE_POSITIONS", common_data_mortar_contact));
350 
351  fe_rhs_mortar_contact->getOpPtrVector().push_back(new OpGetNormalMasterALE(
352  "MESH_NODE_POSITIONS", common_data_mortar_contact));
353 
354  fe_rhs_mortar_contact->getOpPtrVector().push_back(
355  new OpGetPositionAtGaussPtsMaster(field_name,
356  common_data_mortar_contact));
357 
358  fe_rhs_mortar_contact->getOpPtrVector().push_back(
359  new OpGetPositionAtGaussPtsSlave(field_name, common_data_mortar_contact));
360 
361  if (is_eigen_pos_field) {
362  fe_rhs_mortar_contact->getOpPtrVector().push_back(
363  new OpGetDeformationFieldForDisplAtGaussPtsMaster(
364  eigen_pos_field_name, common_data_mortar_contact));
365 
366  fe_rhs_mortar_contact->getOpPtrVector().push_back(
367  new OpGetMatPosForDisplAtGaussPtsMaster("MESH_NODE_POSITIONS",
368  common_data_mortar_contact));
369 
370  fe_rhs_mortar_contact->getOpPtrVector().push_back(
371  new OpGetMatPosForDisplAtGaussPtsSlave("MESH_NODE_POSITIONS",
372  common_data_mortar_contact));
373 
374  fe_rhs_mortar_contact->getOpPtrVector().push_back(
375  new OpGetDeformationFieldForDisplAtGaussPtsSlave(
376  eigen_pos_field_name, common_data_mortar_contact));
377  }
378 
379  if (use_reference_coordinates) {
380  fe_rhs_mortar_contact->getOpPtrVector().push_back(
381  new OpGetDeformationFieldForDisplAtGaussPtsMaster(
382  "MESH_NODE_POSITIONS", common_data_mortar_contact));
383 
384  fe_rhs_mortar_contact->getOpPtrVector().push_back(
385  new OpGetDeformationFieldForDisplAtGaussPtsSlave(
386  "MESH_NODE_POSITIONS", common_data_mortar_contact));
387  }
388 
389  fe_rhs_mortar_contact->getOpPtrVector().push_back(
390  new OpGetGapSlave(field_name, common_data_mortar_contact));
391 
392  fe_rhs_mortar_contact->getOpPtrVector().push_back(
393  new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
394  common_data_mortar_contact));
395 
396  fe_rhs_mortar_contact->getOpPtrVector().push_back(new OpGetGaussPtsState(
397  lagrange_field_name, common_data_mortar_contact, cnValue, is_alm));
398 
399  if (!is_alm) {
400  fe_rhs_mortar_contact->getOpPtrVector().push_back(
401  new OpCalContactTractionOnSlave(field_name,
402  common_data_mortar_contact));
403 
404  fe_rhs_mortar_contact->getOpPtrVector().push_back(new OpCalIntCompFunSlave(
405  lagrange_field_name, common_data_mortar_contact, cnValuePtr));
406  } else {
407  fe_rhs_mortar_contact->getOpPtrVector().push_back(
408  new OpGetAugmentedLambdaSlave(field_name, common_data_mortar_contact,
409  cnValue));
410 
411  fe_rhs_mortar_contact->getOpPtrVector().push_back(
412  new OpCalAugmentedTractionRhsSlave(field_name,
413  common_data_mortar_contact));
414 
415  fe_rhs_mortar_contact->getOpPtrVector().push_back(
416  new OpGapConstraintAugmentedRhs(lagrange_field_name,
417  common_data_mortar_contact, cnValue));
418  }
419 
420  MoFEMFunctionReturn(0);
421 }
422 
423 MoFEMErrorCode MortarContactProblem::setMasterForceOperatorsRhs(
424  boost::shared_ptr<MortarContactElement> fe_rhs_mortar_contact,
425  boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
426  string field_name, string lagrange_field_name, bool is_alm,
427  bool is_eigen_pos_field, string eigen_pos_field_name,
428  bool use_reference_coordinates) {
429  MoFEMFunctionBegin;
430 
431  fe_rhs_mortar_contact->getOpPtrVector().push_back(new OpGetNormalSlaveALE(
432  "MESH_NODE_POSITIONS", common_data_mortar_contact));
433 
434  fe_rhs_mortar_contact->getOpPtrVector().push_back(new OpGetNormalMasterALE(
435  "MESH_NODE_POSITIONS", common_data_mortar_contact));
436 
437  fe_rhs_mortar_contact->getOpPtrVector().push_back(
438  new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
439  common_data_mortar_contact));
440 
441  // fe_rhs_mortar_contact->getOpPtrVector().push_back(
442  // new OpGetPositionAtGaussPtsMaster(field_name,
443  // common_data_mortar_contact));
444 
445  // fe_rhs_mortar_contact->getOpPtrVector().push_back(
446  // new OpGetPositionAtGaussPtsSlave(field_name,
447  // common_data_mortar_contact));
448 
449  if (!is_alm) {
450  fe_rhs_mortar_contact->getOpPtrVector().push_back(
451  new OpCalContactTractionOnMaster(field_name,
452  common_data_mortar_contact));
453  } else {
454 
455  fe_rhs_mortar_contact->getOpPtrVector().push_back(
456  new OpGetPositionAtGaussPtsMaster(field_name,
457  common_data_mortar_contact));
458 
459  fe_rhs_mortar_contact->getOpPtrVector().push_back(
460  new OpGetPositionAtGaussPtsSlave(field_name,
461  common_data_mortar_contact));
462 
463  if (is_eigen_pos_field) {
464  fe_rhs_mortar_contact->getOpPtrVector().push_back(
465  new OpGetDeformationFieldForDisplAtGaussPtsMaster(
466  eigen_pos_field_name, common_data_mortar_contact));
467 
468  fe_rhs_mortar_contact->getOpPtrVector().push_back(
469  new OpGetMatPosForDisplAtGaussPtsMaster("MESH_NODE_POSITIONS",
470  common_data_mortar_contact));
471 
472  fe_rhs_mortar_contact->getOpPtrVector().push_back(
473  new OpGetMatPosForDisplAtGaussPtsSlave("MESH_NODE_POSITIONS",
474  common_data_mortar_contact));
475 
476  fe_rhs_mortar_contact->getOpPtrVector().push_back(
477  new OpGetDeformationFieldForDisplAtGaussPtsSlave(
478  eigen_pos_field_name, common_data_mortar_contact));
479  }
480 
481  if (use_reference_coordinates) {
482  fe_rhs_mortar_contact->getOpPtrVector().push_back(
483  new OpGetDeformationFieldForDisplAtGaussPtsMaster(
484  "MESH_NODE_POSITIONS", common_data_mortar_contact));
485 
486  fe_rhs_mortar_contact->getOpPtrVector().push_back(
487  new OpGetDeformationFieldForDisplAtGaussPtsSlave(
488  "MESH_NODE_POSITIONS", common_data_mortar_contact));
489  }
490 
491  fe_rhs_mortar_contact->getOpPtrVector().push_back(
492  new OpGetGapSlave(field_name, common_data_mortar_contact));
493 
494  fe_rhs_mortar_contact->getOpPtrVector().push_back(
495  new OpGetAugmentedLambdaSlave(field_name, common_data_mortar_contact,
496  cnValue));
497 
498  fe_rhs_mortar_contact->getOpPtrVector().push_back(
499  new OpCalAugmentedTractionRhsMaster(field_name,
500  common_data_mortar_contact));
501  }
502 
503  MoFEMFunctionReturn(0);
504 }
505 
506 MoFEMErrorCode MortarContactProblem::setContactOperatorsLhs(
507  boost::shared_ptr<MortarContactElement> fe_lhs_mortar_contact,
508  boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
509  string field_name, string lagrange_field_name, bool is_alm,
510  bool is_eigen_pos_field, string eigen_pos_field_name,
511  bool use_reference_coordinates) {
512  MoFEMFunctionBegin;
513 
514  fe_lhs_mortar_contact->getOpPtrVector().push_back(
515  new OpGetPositionAtGaussPtsMaster(field_name,
516  common_data_mortar_contact));
517 
518  fe_lhs_mortar_contact->getOpPtrVector().push_back(new OpGetNormalSlaveALE(
519  "MESH_NODE_POSITIONS", common_data_mortar_contact));
520 
521  fe_lhs_mortar_contact->getOpPtrVector().push_back(new OpGetNormalMasterALE(
522  "MESH_NODE_POSITIONS", common_data_mortar_contact));
523 
524  fe_lhs_mortar_contact->getOpPtrVector().push_back(
525  new OpGetPositionAtGaussPtsSlave(field_name, common_data_mortar_contact));
526 
527  if (is_eigen_pos_field) {
528  fe_lhs_mortar_contact->getOpPtrVector().push_back(
529  new OpGetDeformationFieldForDisplAtGaussPtsMaster(
530  eigen_pos_field_name, common_data_mortar_contact));
531 
532  fe_lhs_mortar_contact->getOpPtrVector().push_back(
533  new OpGetMatPosForDisplAtGaussPtsMaster("MESH_NODE_POSITIONS",
534  common_data_mortar_contact));
535 
536  fe_lhs_mortar_contact->getOpPtrVector().push_back(
537  new OpGetMatPosForDisplAtGaussPtsSlave("MESH_NODE_POSITIONS",
538  common_data_mortar_contact));
539 
540  fe_lhs_mortar_contact->getOpPtrVector().push_back(
541  new OpGetDeformationFieldForDisplAtGaussPtsSlave(
542  eigen_pos_field_name, common_data_mortar_contact));
543  }
544 
545  if (use_reference_coordinates) {
546  fe_lhs_mortar_contact->getOpPtrVector().push_back(
547  new OpGetDeformationFieldForDisplAtGaussPtsMaster(
548  "MESH_NODE_POSITIONS", common_data_mortar_contact));
549 
550  fe_lhs_mortar_contact->getOpPtrVector().push_back(
551  new OpGetDeformationFieldForDisplAtGaussPtsSlave(
552  "MESH_NODE_POSITIONS", common_data_mortar_contact));
553  }
554 
555  fe_lhs_mortar_contact->getOpPtrVector().push_back(
556  new OpGetGapSlave(field_name, common_data_mortar_contact));
557 
558  fe_lhs_mortar_contact->getOpPtrVector().push_back(
559  new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
560  common_data_mortar_contact));
561  if (!is_alm) {
562  fe_lhs_mortar_contact->getOpPtrVector().push_back(
563  new OpCalContactTractionOverLambdaSlaveSlave(
564  field_name, lagrange_field_name, common_data_mortar_contact));
565 
566  fe_lhs_mortar_contact->getOpPtrVector().push_back(
567  new OpCalDerIntCompFunOverLambdaSlaveSlave(
568  lagrange_field_name, common_data_mortar_contact, cnValuePtr));
569 
570  fe_lhs_mortar_contact->getOpPtrVector().push_back(
571  new OpCalDerIntCompFunOverSpatPosSlaveMaster(
572  lagrange_field_name, field_name, common_data_mortar_contact,
573  cnValuePtr));
574 
575  fe_lhs_mortar_contact->getOpPtrVector().push_back(
576  new OpCalDerIntCompFunOverSpatPosSlaveSlave(
577  lagrange_field_name, field_name, common_data_mortar_contact,
578  cnValuePtr));
579  } else {
580  fe_lhs_mortar_contact->getOpPtrVector().push_back(
581  new OpGetAugmentedLambdaSlave(field_name, common_data_mortar_contact,
582  cnValue));
583 
584  fe_lhs_mortar_contact->getOpPtrVector().push_back(
585  new OpCalContactAugmentedTractionOverLambdaSlaveSlave(
586  field_name, lagrange_field_name, common_data_mortar_contact));
587 
588  fe_lhs_mortar_contact->getOpPtrVector().push_back(
589  new OpCalContactAugmentedTractionOverSpatialSlaveSlave(
590  field_name, field_name, cnValue, common_data_mortar_contact));
591 
592  fe_lhs_mortar_contact->getOpPtrVector().push_back(
593  new OpCalContactAugmentedTractionOverSpatialSlaveMaster(
594  field_name, field_name, cnValue, common_data_mortar_contact));
595 
596  fe_lhs_mortar_contact->getOpPtrVector().push_back(
597  new OpGapConstraintAugmentedOverLambda(
598  lagrange_field_name, common_data_mortar_contact, cnValue));
599 
600  fe_lhs_mortar_contact->getOpPtrVector().push_back(
601  new OpGapConstraintAugmentedOverSpatialMaster(
602  field_name, lagrange_field_name, common_data_mortar_contact,
603  cnValue));
604 
605  fe_lhs_mortar_contact->getOpPtrVector().push_back(
606  new OpGapConstraintAugmentedOverSpatialSlave(
607  field_name, lagrange_field_name, common_data_mortar_contact,
608  cnValue));
609  }
610  MoFEMFunctionReturn(0);
611 }
612 
613 MoFEMErrorCode MortarContactProblem::setContactOperatorsLhs(
614  boost::shared_ptr<MortarConvectMasterContactElement> fe_lhs_mortar_contact,
615  boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
616  string field_name, string lagrange_field_name, bool is_alm,
617  bool is_eigen_pos_field, string eigen_pos_field_name,
618  bool use_reference_coordinates) {
619  MoFEMFunctionBegin;
620 
621  CHKERR MortarContactProblem::setContactOperatorsLhs(
622  boost::dynamic_pointer_cast<MortarContactElement>(fe_lhs_mortar_contact),
623  common_data_mortar_contact, field_name, lagrange_field_name, is_alm,
624  is_eigen_pos_field, eigen_pos_field_name, use_reference_coordinates);
625 
626  // CHKERR setContactOperatorsLhs(
627  // boost::dynamic_pointer_cast<MortarContactElement>(
628  // fe_lhs_mortar_contact),
629  // common_data_mortar_contact, field_name, lagrange_field_name);
630 
631  fe_lhs_mortar_contact->getOpPtrVector().push_back(
632  new OpCalculateGradPositionXi(field_name, common_data_mortar_contact));
633 
634  fe_lhs_mortar_contact->getOpPtrVector().push_back(
635  new OpLhsConvectIntegrationPtsConstrainMasterGap(
636  lagrange_field_name, field_name, common_data_mortar_contact,
637  cnValuePtr, ContactOp::FACESLAVESLAVE,
638  fe_lhs_mortar_contact->getConvectPtr()->getDiffKsiSpatialSlave()));
639 
640  fe_lhs_mortar_contact->getOpPtrVector().push_back(
641  new OpLhsConvectIntegrationPtsConstrainMasterGap(
642  lagrange_field_name, field_name, common_data_mortar_contact,
643  cnValuePtr, ContactOp::FACESLAVEMASTER,
644  fe_lhs_mortar_contact->getConvectPtr()->getDiffKsiSpatialMaster()));
645 
646  MoFEMFunctionReturn(0);
647 }
648 
649 MoFEMErrorCode MortarContactProblem::setMasterForceOperatorsLhs(
650  boost::shared_ptr<MortarContactElement> fe_lhs_mortar_contact,
651  boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
652  string field_name, string lagrange_field_name, bool is_alm,
653  bool is_eigen_pos_field, string eigen_pos_field_name,
654  bool use_reference_coordinates) {
655  MoFEMFunctionBegin;
656 
657  fe_lhs_mortar_contact->getOpPtrVector().push_back(new OpGetNormalSlaveALE(
658  "MESH_NODE_POSITIONS", common_data_mortar_contact));
659 
660  fe_lhs_mortar_contact->getOpPtrVector().push_back(new OpGetNormalMasterALE(
661  "MESH_NODE_POSITIONS", common_data_mortar_contact));
662 
663  fe_lhs_mortar_contact->getOpPtrVector().push_back(
664  new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
665  common_data_mortar_contact));
666  if (!is_alm) {
667  fe_lhs_mortar_contact->getOpPtrVector().push_back(
668  new OpCalContactTractionOverLambdaMasterSlave(
669  field_name, lagrange_field_name, common_data_mortar_contact));
670  } else {
671 
672  fe_lhs_mortar_contact->getOpPtrVector().push_back(
673  new OpGetPositionAtGaussPtsMaster(field_name,
674  common_data_mortar_contact));
675 
676  fe_lhs_mortar_contact->getOpPtrVector().push_back(
677  new OpGetPositionAtGaussPtsSlave(field_name,
678  common_data_mortar_contact));
679 
680  if (is_eigen_pos_field) {
681  fe_lhs_mortar_contact->getOpPtrVector().push_back(
682  new OpGetDeformationFieldForDisplAtGaussPtsMaster(
683  eigen_pos_field_name, common_data_mortar_contact));
684 
685  fe_lhs_mortar_contact->getOpPtrVector().push_back(
686  new OpGetMatPosForDisplAtGaussPtsMaster("MESH_NODE_POSITIONS",
687  common_data_mortar_contact));
688 
689  fe_lhs_mortar_contact->getOpPtrVector().push_back(
690  new OpGetMatPosForDisplAtGaussPtsSlave("MESH_NODE_POSITIONS",
691  common_data_mortar_contact));
692 
693  fe_lhs_mortar_contact->getOpPtrVector().push_back(
694  new OpGetDeformationFieldForDisplAtGaussPtsSlave(
695  eigen_pos_field_name, common_data_mortar_contact));
696  }
697 
698  if (use_reference_coordinates) {
699  fe_lhs_mortar_contact->getOpPtrVector().push_back(
700  new OpGetDeformationFieldForDisplAtGaussPtsMaster(
701  "MESH_NODE_POSITIONS", common_data_mortar_contact));
702 
703  fe_lhs_mortar_contact->getOpPtrVector().push_back(
704  new OpGetDeformationFieldForDisplAtGaussPtsSlave(
705  "MESH_NODE_POSITIONS", common_data_mortar_contact));
706  }
707 
708  fe_lhs_mortar_contact->getOpPtrVector().push_back(
709  new OpGetGapSlave(field_name, common_data_mortar_contact));
710 
711  fe_lhs_mortar_contact->getOpPtrVector().push_back(
712  new OpGetAugmentedLambdaSlave(field_name, common_data_mortar_contact,
713  cnValue));
714 
715  fe_lhs_mortar_contact->getOpPtrVector().push_back(
716  new OpCalContactAugmentedTractionOverLambdaMasterSlave(
717  field_name, lagrange_field_name, common_data_mortar_contact));
718 
719  fe_lhs_mortar_contact->getOpPtrVector().push_back(
720  new OpCalContactAugmentedTractionOverSpatialMasterMaster(
721  field_name, field_name, cnValue, common_data_mortar_contact));
722 
723  fe_lhs_mortar_contact->getOpPtrVector().push_back(
724  new OpCalContactAugmentedTractionOverSpatialMasterSlave(
725  field_name, field_name, cnValue, common_data_mortar_contact));
726  }
727 
728  MoFEMFunctionReturn(0);
729 }
730 
731 MoFEMErrorCode MortarContactProblem::setMasterForceOperatorsLhs(
732  boost::shared_ptr<MortarConvectSlaveContactElement> fe_lhs_mortar_contact,
733  boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
734  string field_name, string lagrange_field_name, bool is_alm,
735  bool is_eigen_pos_field, string eigen_pos_field_name,
736  bool use_reference_coordinates) {
737  MoFEMFunctionBegin;
738 
739  CHKERR MortarContactProblem::setMasterForceOperatorsLhs(
740  boost::dynamic_pointer_cast<MortarContactElement>(fe_lhs_mortar_contact),
741  common_data_mortar_contact, field_name, lagrange_field_name, is_alm,
742  is_eigen_pos_field, eigen_pos_field_name, use_reference_coordinates);
743 
744  fe_lhs_mortar_contact->getOpPtrVector().push_back(new OpCalculateGradLambdaXi(
745  lagrange_field_name, common_data_mortar_contact));
746 
747  fe_lhs_mortar_contact->getOpPtrVector().push_back(
748  new OpLhsConvectIntegrationPtsContactTraction(
749  field_name, field_name, common_data_mortar_contact,
750  ContactOp::FACEMASTERSLAVE,
751  fe_lhs_mortar_contact->getConvectPtr()->getDiffKsiSpatialSlave()));
752 
753  fe_lhs_mortar_contact->getOpPtrVector().push_back(
754  new OpLhsConvectIntegrationPtsContactTraction(
755  field_name, field_name, common_data_mortar_contact,
756  ContactOp::FACEMASTERMASTER,
757  fe_lhs_mortar_contact->getConvectPtr()->getDiffKsiSpatialMaster()));
758  MoFEMFunctionReturn(0);
759 }
760 
761 MoFEMErrorCode MortarContactProblem::setContactOperatorsForPostProc(
762  boost::shared_ptr<MortarContactElement> fe_post_proc_mortar_contact,
763  boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
764  MoFEM::Interface &m_field, string field_name, string lagrange_field_name,
765  moab::Interface &moab_out, bool alm_flag, bool is_eigen_pos_field,
766  string eigen_pos_field_name, bool is_displacements,
767  Range post_proc_surface, double post_proc_gap_tol) {
768  MoFEMFunctionBegin;
769 
770  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
771  new OpGetNormalMasterALE("MESH_NODE_POSITIONS",
772  common_data_mortar_contact));
773 
774  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
775  new OpGetNormalSlaveALE("MESH_NODE_POSITIONS",
776  common_data_mortar_contact));
777 
778  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
779  new OpGetPositionAtGaussPtsMaster(field_name,
780  common_data_mortar_contact));
781 
782  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
783  new OpGetPositionAtGaussPtsSlave(field_name, common_data_mortar_contact));
784 
785  if (is_eigen_pos_field) {
786  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
787  new OpGetDeformationFieldForDisplAtGaussPtsMaster(
788  eigen_pos_field_name, common_data_mortar_contact));
789 
790  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
791  new OpGetMatPosForDisplAtGaussPtsMaster("MESH_NODE_POSITIONS",
792  common_data_mortar_contact));
793 
794  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
795  new OpGetMatPosForDisplAtGaussPtsSlave("MESH_NODE_POSITIONS",
796  common_data_mortar_contact));
797 
798  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
799  new OpGetDeformationFieldForDisplAtGaussPtsSlave(
800  eigen_pos_field_name, common_data_mortar_contact));
801  }
802 
803  if (is_displacements && !is_eigen_pos_field) {
804  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
805  new OpGetDeformationFieldForDisplAtGaussPtsMaster(
806  "MESH_NODE_POSITIONS", common_data_mortar_contact));
807 
808  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
809  new OpGetDeformationFieldForDisplAtGaussPtsSlave(
810  "MESH_NODE_POSITIONS", common_data_mortar_contact));
811  }
812  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
813  new OpGetGapSlave(field_name, common_data_mortar_contact));
814 
815  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
816  new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
817  common_data_mortar_contact));
818 
819  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
820  new OpLagGapProdGaussPtsSlave(lagrange_field_name,
821  common_data_mortar_contact));
822 
823  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
824  new OpGetGaussPtsState(lagrange_field_name, common_data_mortar_contact,
825  cnValue, alm_flag));
826 
827  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
828  new OpMakeVtkSlave(m_field, field_name, common_data_mortar_contact,
829  moab_out, NO_TAG, post_proc_surface));
830 
831  fe_post_proc_mortar_contact->getOpPtrVector().push_back(new OpGetContactArea(
832  lagrange_field_name, common_data_mortar_contact, cnValue, alm_flag,
833  post_proc_surface, post_proc_gap_tol));
834 
835  MoFEMFunctionReturn(0);
836 }
837 
838 // setup operators for calculation of active set
839 MoFEMErrorCode MortarContactProblem::setContactOperatorsRhsALEMaterial(
840  boost::shared_ptr<MortarContactElement> fe_rhs_mortar_contact_ale,
841  boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
842  const string field_name, const string mesh_node_field_name,
843  const string lagrange_field_name, const string side_fe_name) {
844  MoFEMFunctionBegin;
845 
846  boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnContactPrismSide>
847  fe_mat_side_rhs_master = boost::make_shared<
848  VolumeElementForcesAndSourcesCoreOnContactPrismSide>(mField);
849 
850  fe_mat_side_rhs_master->getOpPtrVector().push_back(
851  new OpCalculateVectorFieldGradient<3, 3>(
852  mesh_node_field_name, common_data_mortar_contact->HMat));
853  fe_mat_side_rhs_master->getOpPtrVector().push_back(
854  new OpCalculateVectorFieldGradient<3, 3>(
855  field_name, common_data_mortar_contact->hMat));
856 
857  boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnContactPrismSide>
858  fe_mat_side_rhs_slave = boost::make_shared<
859  VolumeElementForcesAndSourcesCoreOnContactPrismSide>(mField);
860 
861  fe_mat_side_rhs_slave->getOpPtrVector().push_back(
862  new OpCalculateVectorFieldGradient<3, 3>(
863  mesh_node_field_name, common_data_mortar_contact->HMat));
864  fe_mat_side_rhs_slave->getOpPtrVector().push_back(
865  new OpCalculateVectorFieldGradient<3, 3>(
866  field_name, common_data_mortar_contact->hMat));
867 
868  fe_rhs_mortar_contact_ale->getOpPtrVector().push_back(new OpGetNormalSlaveALE(
869  "MESH_NODE_POSITIONS", common_data_mortar_contact));
870 
871  fe_rhs_mortar_contact_ale->getOpPtrVector().push_back(
872  new OpGetNormalMasterALE("MESH_NODE_POSITIONS",
873  common_data_mortar_contact));
874 
875  fe_rhs_mortar_contact_ale->getOpPtrVector().push_back(
876  new OpGetPositionAtGaussPtsMaster(field_name,
877  common_data_mortar_contact));
878 
879  fe_rhs_mortar_contact_ale->getOpPtrVector().push_back(
880  new OpGetPositionAtGaussPtsSlave(field_name, common_data_mortar_contact));
881 
882  fe_rhs_mortar_contact_ale->getOpPtrVector().push_back(
883  new OpGetGapSlave(field_name, common_data_mortar_contact));
884 
885  fe_rhs_mortar_contact_ale->getOpPtrVector().push_back(
886  new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
887  common_data_mortar_contact));
888 
889  fe_mat_side_rhs_master->getOpPtrVector().push_back(new OpCalculateDeformation(
890  mesh_node_field_name, common_data_mortar_contact, false));
891 
892  fe_rhs_mortar_contact_ale->getOpPtrVector().push_back(
893  new OpLoopForSideOfContactPrism(mesh_node_field_name,
894  fe_mat_side_rhs_master, side_fe_name,
895  ContactOp::FACEMASTER));
896 
897  fe_rhs_mortar_contact_ale->getOpPtrVector().push_back(
898  new OpCalMatForcesALEMaster(mesh_node_field_name,
899  common_data_mortar_contact));
900 
901  fe_mat_side_rhs_slave->getOpPtrVector().push_back(new OpCalculateDeformation(
902  mesh_node_field_name, common_data_mortar_contact, false));
903 
904  fe_rhs_mortar_contact_ale->getOpPtrVector().push_back(
905  new OpLoopForSideOfContactPrism(mesh_node_field_name,
906  fe_mat_side_rhs_slave, side_fe_name,
907  ContactOp::FACESLAVE));
908 
909  fe_rhs_mortar_contact_ale->getOpPtrVector().push_back(
910  new OpCalMatForcesALESlave(mesh_node_field_name,
911  common_data_mortar_contact));
912  MoFEMFunctionReturn(0);
913 }
914 
915 MoFEMErrorCode MortarContactProblem::setContactOperatorsLhsALEMaterial(
916  boost::shared_ptr<MortarContactElement> fe_lhs_mortar_contact_ale,
917  boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
918  const string field_name, const string mesh_node_field_name,
919  const string lagrange_field_name, const string side_fe_name) {
920  MoFEMFunctionBegin;
921 
922  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(new OpGetNormalSlaveALE(
923  mesh_node_field_name, common_data_mortar_contact));
924 
925  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
926  new OpGetNormalMasterALE(mesh_node_field_name,
927  common_data_mortar_contact));
928 
929  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
930  new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
931  common_data_mortar_contact));
932 
933  boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnContactPrismSide>
934  feMatSideLhs_dx = boost::make_shared<
935  VolumeElementForcesAndSourcesCoreOnContactPrismSide>(mField);
936 
937  feMatSideLhs_dx->getOpPtrVector().push_back(
938  new OpCalculateVectorFieldGradient<3, 3>(
939  mesh_node_field_name, common_data_mortar_contact->HMat));
940 
941  feMatSideLhs_dx->getOpPtrVector().push_back(
942  new OpCalculateVectorFieldGradient<3, 3>(
943  field_name, common_data_mortar_contact->hMat));
944 
945  // // Master derivative over spatial
946  feMatSideLhs_dx->getOpPtrVector().push_back(new OpCalculateDeformation(
947  mesh_node_field_name, common_data_mortar_contact, false));
948 
949  feMatSideLhs_dx->getOpPtrVector().push_back(
950  new OpContactMaterialVolOnSideLhs_dX_dX(
951  mesh_node_field_name, mesh_node_field_name,
952  common_data_mortar_contact, true));
953 
954  feMatSideLhs_dx->getOpPtrVector().push_back(
955  new OpContactMaterialVolOnSideLhs_dX_dx(
956  mesh_node_field_name, field_name, common_data_mortar_contact, true));
957 
958  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
959  new OpLoopForSideOfContactPrism(mesh_node_field_name, feMatSideLhs_dx,
960  side_fe_name, ContactOp::FACEMASTER));
961 
962  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
963  new OpContactMaterialMasterOnFaceLhs_dX_dX(
964  mesh_node_field_name, mesh_node_field_name,
965  common_data_mortar_contact, POSITION_RANK, POSITION_RANK));
966 
967  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
968  new OpContactMaterialMasterSlaveLhs_dX_dLagmult(
969  mesh_node_field_name, lagrange_field_name, common_data_mortar_contact,
970  POSITION_RANK, LAGRANGE_RANK));
971 
972  boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnContactPrismSide>
973  feMatSideLhsSlave_dx = boost::make_shared<
974  VolumeElementForcesAndSourcesCoreOnContactPrismSide>(mField);
975 
976  feMatSideLhsSlave_dx->getOpPtrVector().push_back(
977  new OpCalculateVectorFieldGradient<3, 3>(
978  mesh_node_field_name, common_data_mortar_contact->HMat));
979 
980  feMatSideLhsSlave_dx->getOpPtrVector().push_back(
981  new OpCalculateVectorFieldGradient<3, 3>(
982  field_name, common_data_mortar_contact->hMat));
983 
984  feMatSideLhsSlave_dx->getOpPtrVector().push_back(new OpCalculateDeformation(
985  mesh_node_field_name, common_data_mortar_contact, false));
986 
987  feMatSideLhsSlave_dx->getOpPtrVector().push_back(
988  new OpContactMaterialVolOnSideLhs_dX_dX(
989  mesh_node_field_name, mesh_node_field_name,
990  common_data_mortar_contact, false));
991 
992  feMatSideLhsSlave_dx->getOpPtrVector().push_back(
993  new OpContactMaterialVolOnSideLhs_dX_dx(
994  mesh_node_field_name, field_name, common_data_mortar_contact, false));
995 
996  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
997  new OpLoopForSideOfContactPrism(mesh_node_field_name,
998  feMatSideLhsSlave_dx, side_fe_name,
999  ContactOp::FACESLAVE));
1000 
1001  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
1002  new OpContactMaterialSlaveOnFaceLhs_dX_dX(
1003  mesh_node_field_name, mesh_node_field_name,
1004  common_data_mortar_contact, POSITION_RANK, POSITION_RANK));
1005 
1006  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
1007  new OpContactMaterialSlaveSlaveLhs_dX_dLagmult(
1008  mesh_node_field_name, lagrange_field_name, common_data_mortar_contact,
1009  POSITION_RANK, LAGRANGE_RANK));
1010 
1011  MoFEMFunctionReturn(0);
1012 }
1013 
1014 MoFEMErrorCode MortarContactProblem::setContactOperatorsLhsALE(
1015  boost::shared_ptr<MortarContactElement> fe_lhs_mortar_contact_ale,
1016  boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
1017  const string field_name, const string mesh_node_field_name,
1018  const string lagrange_field_name) {
1019  MoFEMFunctionBegin;
1020 
1021  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(new OpGetNormalSlaveALE(
1022  mesh_node_field_name, common_data_mortar_contact));
1023 
1024  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
1025  new OpGetNormalMasterALE(mesh_node_field_name,
1026  common_data_mortar_contact));
1027 
1028  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
1029  new OpGetPositionAtGaussPtsMaster(field_name,
1030  common_data_mortar_contact));
1031 
1032  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
1033  new OpGetPositionAtGaussPtsSlave(field_name, common_data_mortar_contact));
1034 
1035  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
1036  new OpGetGapSlave(field_name, common_data_mortar_contact));
1037 
1038  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
1039  new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
1040  common_data_mortar_contact));
1041 
1042  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
1043  new OpContactTractionSlaveSlave_dX(field_name, mesh_node_field_name,
1044  common_data_mortar_contact,
1045  POSITION_RANK, POSITION_RANK));
1046 
1047  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
1048  new OpContactTractionMasterSlave_dX(field_name, mesh_node_field_name,
1049  common_data_mortar_contact,
1050  POSITION_RANK, POSITION_RANK));
1051 
1052  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
1053  new OpContactTractionMasterMaster_dX(field_name, mesh_node_field_name,
1054  common_data_mortar_contact,
1055  POSITION_RANK, POSITION_RANK));
1056 
1057  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
1058  new OpCalDerIntCompFunSlaveSlave_dX(
1059  lagrange_field_name, mesh_node_field_name, cnValuePtr,
1060  common_data_mortar_contact, LAGRANGE_RANK, POSITION_RANK));
1061  MoFEMFunctionReturn(0);
1062 }
SimpleContactProblem::OpContactMaterialSlaveSlaveLhs_dX_dLagmult
LHS-operator for the contact element (material configuration)
Definition: SimpleContact.hpp:2944
SimpleContactProblem::OpCalContactAugmentedTractionOverLambdaSlaveSlave
LHS-operator for the simple contact element with Augmented Lagrangian Method.
Definition: SimpleContact.hpp:1642
SimpleContactProblem::OpCalculateDeformation
Operator for computing deformation gradients in side volumes.
Definition: SimpleContact.hpp:2456
SimpleContactProblem::OpGapConstraintAugmentedOverSpatialSlave
LHS-operator for the simple contact element.
Definition: SimpleContact.hpp:2116
MortarContactProblem::setMasterForceOperatorsRhs
MoFEMErrorCode setMasterForceOperatorsRhs(boost::shared_ptr< MortarContactElement > fe_rhs_mortar_contact, boost::shared_ptr< CommonDataMortarContact > common_data_mortar_contact, string field_name, string lagrange_field_name, bool is_alm=false, bool is_eigen_pos_field=false, string eigen_pos_field_name="EIGEN_SPATIAL_POSITIONS", bool use_reference_coordinates=false)
Definition: MortarContactProblem.cpp:423
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
CHKERRQ_MOAB
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:467
SimpleContactProblem::OpCalDerIntCompFunSlaveSlave_dX
LHS-operator for the simple contact element.
Definition: SimpleContact.hpp:3236
NOBASE
@ NOBASE
Definition: definitions.h:59
QUAD_2D_TABLE
static QUAD *const QUAD_2D_TABLE[]
Definition: quad.h:175
SimpleContactProblem::OpGetMatPosForDisplAtGaussPtsMaster
Operator for the simple contact element.
Definition: SimpleContact.hpp:773
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
SimpleContactProblem::OpContactMaterialVolOnSideLhs_dX_dX
LHS-operator for the contact element (material configuration)
Definition: SimpleContact.hpp:3397
SimpleContactProblem::OpGetNormalSlaveALE
Computes, for material configuration, normal to slave face that is common to all gauss points.
Definition: SimpleContact.hpp:2645
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
SimpleContactProblem::OpCalContactAugmentedTractionOverSpatialSlaveMaster
LHS-operator for the simple contact element with Augmented Lagrangian Method.
Definition: SimpleContact.hpp:1917
MortarContactProblem::setContactOperatorsLhsALEMaterial
MoFEMErrorCode setContactOperatorsLhsALEMaterial(boost::shared_ptr< MortarContactElement > fe_lhs_simple_contact_ale, boost::shared_ptr< CommonDataMortarContact > common_data_simple_contact, const string field_name, const string mesh_node_field_name, const string lagrange_field_name, const string side_fe_name)
Function for the mortar contact element that sets the user data LHS-operators.
Definition: MortarContactProblem.cpp:915
ContactSearchKdTree
Definition: ContactSearchKdTree.hpp:18
SimpleContactProblem::OpCalContactTractionOnMaster
RHS-operator for the simple contact element.
Definition: SimpleContact.hpp:962
SimpleContactProblem::OpContactTractionSlaveSlave_dX
LHS-operator for the simple contact element.
Definition: SimpleContact.hpp:3023
MortarContactProblem::MortarContactElement::area2D
double area2D(double *a, double *b, double *c)
Definition: MortarContactProblem.hpp:50
SimpleContactProblem::OpCalDerIntCompFunOverSpatPosSlaveMaster
LHS-operator for the simple contact element.
Definition: SimpleContact.hpp:1448
SimpleContactProblem::OpContactMaterialVolOnSideLhs_dX_dx
LHS-operator (material configuration) on the side volume of either master or slave side.
Definition: SimpleContact.hpp:3341
FTensor::Tensor2< double *, 3, 3 >
order
constexpr int order
Definition: dg_projection.cpp:18
SimpleContactProblem::OpCalMatForcesALESlave
RHS - operator for the contact element (material configuration) Integrates virtual work of contact tr...
Definition: SimpleContact.hpp:2581
SimpleContactProblem::OpCalContactAugmentedTractionOverSpatialMasterSlave
LHS-operator for the simple contact element with Augmented Lagrangian Method.
Definition: SimpleContact.hpp:1776
SimpleContactProblem::OpGapConstraintAugmentedOverSpatialMaster
LHS-operator for the simple contact element.
Definition: SimpleContact.hpp:2051
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
SimpleContactProblem::OpLoopForSideOfContactPrism
Trigers operators for side volume element on slave side for evaluation of the RHS contact traction in...
Definition: SimpleContact.hpp:2487
SimpleContactProblem::OpLhsConvectIntegrationPtsContactTraction
Calculate tangent operator for contact force for change of integration point positions,...
Definition: SimpleContact.hpp:2361
ContactSearchKdTree::rotationMatrix
MoFEMErrorCode rotationMatrix(MatrixDouble &m_rot, VectorDouble &v_coords)
Definition: ContactSearchKdTree.hpp:55
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
SimpleContactProblem::OpCalContactTractionOnSlave
RHS-operator for the simple contact element.
Definition: SimpleContact.hpp:1007
SimpleContactProblem::OpGetContactArea
Definition: SimpleContact.hpp:3543
ShapeDiffMBTRI
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:194
SimpleContactProblem::OpLagGapProdGaussPtsSlave
Operator for the simple contact element.
Definition: SimpleContact.hpp:942
SimpleContactProblem::OpGetGapSlave
Operator for the simple contact element.
Definition: SimpleContact.hpp:863
SimpleContactProblem::OpGapConstraintAugmentedOverLambda
LHS-operator for the simple contact element.
Definition: SimpleContact.hpp:1987
SimpleContactProblem::OpGetLagMulAtGaussPtsSlave
Operator for the simple contact element.
Definition: SimpleContact.hpp:897
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MortarContactProblem::setContactOperatorsLhs
MoFEMErrorCode setContactOperatorsLhs(boost::shared_ptr< MortarContactElement > fe_lhs_mortar_contact, boost::shared_ptr< CommonDataMortarContact > common_data_mortar_contact, string field_name, string lagrange_field_name, bool is_alm=false, bool is_eigen_pos_field=false, string eigen_pos_field_name="EIGEN_SPATIAL_POSITIONS", bool use_reference_coordinates=false)
Definition: MortarContactProblem.cpp:506
MortarContactProblem::setContactOperatorsRhs
MoFEMErrorCode setContactOperatorsRhs(boost::shared_ptr< MortarContactElement > fe_rhs_mortar_contact, boost::shared_ptr< CommonDataMortarContact > common_data_mortar_contact, string field_name, string lagrange_field_name, bool is_alm=false, bool is_eigen_pos_field=false, string eigen_pos_field_name="EIGEN_SPATIAL_POSITIONS", bool use_reference_coordinates=false)
Definition: MortarContactProblem.cpp:340
SimpleContactProblem::OpCalContactTractionOverLambdaSlaveSlave
LHS-operator for the simple contact element.
Definition: SimpleContact.hpp:1332
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
SimpleContactProblem::OpGapConstraintAugmentedRhs
RHS-operator for the simple contact element for Augmented Lagrangian Method.
Definition: SimpleContact.hpp:1214
MortarContactProblem::MortarConvectSlaveContactElement::setGaussPts
MoFEMErrorCode setGaussPts(int order)
Definition: MortarContactProblem.cpp:333
MortarContactProblem::MortarContactElement::setGaussPts
virtual MoFEMErrorCode setGaussPts(int order)
Definition: MortarContactProblem.cpp:16
SimpleContactProblem::OpContactMaterialMasterOnFaceLhs_dX_dX
LHS-operator for the contact element (material configuration)
Definition: SimpleContact.hpp:2739
SimpleContactProblem::OpContactTractionMasterSlave_dX
LHS-operator for the simple contact element.
Definition: SimpleContact.hpp:3086
SimpleContactProblem::NO_TAG
@ NO_TAG
Definition: SimpleContact.hpp:2177
ShapeMBTRI_inverse
PetscErrorCode ShapeMBTRI_inverse(double *N, double *diffN, const double *elem_coords, const double *glob_coords, double *loc_coords)
calculate local coordinates of triangle element for given global coordinates in 2D (Assume e....
Definition: fem_tools.c:380
SimpleContactProblem::OpGetPositionAtGaussPtsMaster
Operator for the simple contact element.
Definition: SimpleContact.hpp:755
SimpleContactProblem::OpCalDerIntCompFunOverSpatPosSlaveSlave
LHS-operator for the simple contact element.
Definition: SimpleContact.hpp:1512
SimpleContactProblem::OpCalAugmentedTractionRhsMaster
RHS-operator for the simple contact element.
Definition: SimpleContact.hpp:1051
ContactSearchKdTree::Prism_tag
Definition: ContactSearchKdTree.hpp:44
SimpleContactProblem::OpCalDerIntCompFunOverLambdaSlaveSlave
LHS-operator for the simple contact element.
Definition: SimpleContact.hpp:1386
QUAD_::npoints
int npoints
Definition: quad.h:29
MortarContactProblem::cnValue
double cnValue
Definition: MortarContactProblem.hpp:256
SimpleContactProblem::OpGetPositionAtGaussPtsSlave
Operator for the simple contact element.
Definition: SimpleContact.hpp:809
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
SimpleContactProblem::OpCalContactAugmentedTractionOverSpatialMasterMaster
LHS-operator for the simple contact element with Augmented Lagrangian Method.
Definition: SimpleContact.hpp:1706
MortarContactProblem::setContactOperatorsLhsALE
MoFEMErrorCode setContactOperatorsLhsALE(boost::shared_ptr< MortarContactElement > fe_lhs_simple_contact_ale, boost::shared_ptr< CommonDataMortarContact > common_data_simple_contact, const string field_name, const string mesh_node_field_name, const string lagrange_field_name)
Function for the mortar contact element that sets the user data LHS-operators.
Definition: MortarContactProblem.cpp:1014
SimpleContactProblem::OpCalculateGradLambdaXi
Evaluate gradient of Lagrange multipliers on reference slave surface.
Definition: SimpleContact.hpp:2408
SimpleContactProblem::OpGetMatPosForDisplAtGaussPtsSlave
Operator for the simple contact element.
Definition: SimpleContact.hpp:827
SimpleContactProblem::OpContactMaterialSlaveOnFaceLhs_dX_dX
LHS-operator for the contact element (material configuration)
Definition: SimpleContact.hpp:2802
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
SimpleContactProblem::OpGetDeformationFieldForDisplAtGaussPtsSlave
Operator for the simple contact element.
Definition: SimpleContact.hpp:845
FTensor::Index< 'i', 3 >
MortarContactProblem::MortarContactElement::contactCommondataMultiIndex
boost::shared_ptr< ContactSearchKdTree::ContactCommonData_multiIndex > contactCommondataMultiIndex
Definition: MortarContactProblem.hpp:33
SimpleContactProblem::OpGetAugmentedLambdaSlave
Operator for the simple contact element for Augmented Lagrangian Method.
Definition: SimpleContact.hpp:918
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
SimpleContactProblem::POSITION_RANK
static constexpr int POSITION_RANK
Definition: SimpleContact.hpp:63
Range
SimpleContactProblem::OpCalculateGradPositionXi
Evaluate gradient position on reference master surface.
Definition: SimpleContact.hpp:2390
SimpleContactProblem::LAGRANGE_RANK
static constexpr int LAGRANGE_RANK
Definition: SimpleContact.hpp:62
SimpleContactProblem::OpGetGaussPtsState
Definition: SimpleContact.hpp:3523
SimpleContactProblem::OpCalContactAugmentedTractionOverLambdaMasterSlave
LHS-operator for the simple contact element with Augmented Lagrangian Method.
Definition: SimpleContact.hpp:1577
MortarContactProblem::MortarConvectMasterContactElement::setGaussPts
MoFEMErrorCode setGaussPts(int order)
Definition: MortarContactProblem.cpp:324
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
SimpleContactProblem::OpMakeVtkSlave
Operator for the simple contact element.
Definition: SimpleContact.hpp:2186
SimpleContactProblem::OpContactMaterialMasterSlaveLhs_dX_dLagmult
LHS-operator for the contact element (material configuration)
Definition: SimpleContact.hpp:2864
MortarContactProblem::setContactOperatorsRhsALEMaterial
MoFEMErrorCode setContactOperatorsRhsALEMaterial(boost::shared_ptr< MortarContactElement > fe_rhs_simple_contact_ale, boost::shared_ptr< CommonDataMortarContact > common_data_simple_contact, const string field_name, const string mesh_node_field_name, const string lagrange_field_name, const string side_fe_name)
Function for the mortar contact element that sets the user data post processing operators.
Definition: MortarContactProblem.cpp:839
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
SimpleContactProblem::OpGetDeformationFieldForDisplAtGaussPtsMaster
Operator for the simple contact element.
Definition: SimpleContact.hpp:791
SimpleContactProblem::OpCalAugmentedTractionRhsSlave
RHS-operator for the simple contact element.
Definition: SimpleContact.hpp:1104
MortarContactProblem::cnValuePtr
boost::shared_ptr< double > cnValuePtr
Definition: MortarContactProblem.hpp:257
SimpleContactProblem::OpGetNormalMasterALE
Computes, for material configuration, normal to master face that is common to all gauss points.
Definition: SimpleContact.hpp:2689
SimpleContactProblem::OpCalIntCompFunSlave
RHS-operator for the simple contact element.
Definition: SimpleContact.hpp:1160
MortarContactProblem::MortarContactElement::mField
MoFEM::Interface & mField
Definition: MortarContactProblem.hpp:31
MortarContactProblem::setContactOperatorsForPostProc
MoFEMErrorCode setContactOperatorsForPostProc(boost::shared_ptr< MortarContactElement > fe_post_proc_mortar_contact, boost::shared_ptr< CommonDataMortarContact > common_data_mortar_contact, MoFEM::Interface &m_field, string field_name, string lagrange_field_name, moab::Interface &moab_out, bool alm_flag=false, bool is_eigen_pos_field=false, string eigen_pos_field_name="EIGEN_SPATIAL_POSITIONS", bool is_displacements=false, Range post_proc_surface=Range(), double post_proc_gap_tol=0.0)
Definition: MortarContactProblem.cpp:761
SimpleContactProblem::OpCalContactAugmentedTractionOverSpatialSlaveSlave
LHS-operator for the simple contact element with Augmented Lagrangian Method.
Definition: SimpleContact.hpp:1846
SimpleContactProblem::OpLhsConvectIntegrationPtsConstrainMasterGap
Tangent opeerator for contrains equation for change of spatial positions on master and slave.
Definition: SimpleContact.hpp:2427
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
SimpleContactProblem::OpCalContactTractionOverLambdaMasterSlave
LHS-operator for the simple contact element.
Definition: SimpleContact.hpp:1277
MortarContactProblem::setMasterForceOperatorsLhs
MoFEMErrorCode setMasterForceOperatorsLhs(boost::shared_ptr< MortarContactElement > fe_lhs_mortar_contact, boost::shared_ptr< CommonDataMortarContact > common_data_mortar_contact, string field_name, string lagrange_field_name, bool is_alm=false, bool is_eigen_pos_field=false, string eigen_pos_field_name="EIGEN_SPATIAL_POSITIONS", bool use_reference_coordinates=false)
Definition: MortarContactProblem.cpp:649
OpCalculateVectorFieldGradient
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
SimpleContactProblem::OpContactTractionMasterMaster_dX
LHS-operator for the simple contact element.
Definition: SimpleContact.hpp:3168
ShapeMBTRI
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition: fem_tools.c:182
MortarContactProblem::mField
MoFEM::Interface & mField
Definition: MortarContactProblem.hpp:259
SimpleContactProblem::OpCalMatForcesALEMaster
RHS - operator for the contact element (material configuration) Integrates virtual work of contact tr...
Definition: SimpleContact.hpp:2514
Generated by Doxygen 1.8.17 and hosted at University of Glasgow