v0.14.0
Public Member Functions | Public Attributes | List of all members
MortarContactProblem::MortarContactElement Struct Reference

#include <users_modules/mortar_contact/src/MortarContactProblem.hpp>

Inheritance diagram for MortarContactProblem::MortarContactElement:
[legend]
Collaboration diagram for MortarContactProblem::MortarContactElement:
[legend]

Public Member Functions

 MortarContactElement (MoFEM::Interface &m_field, boost::shared_ptr< ContactSearchKdTree::ContactCommonData_multiIndex > contact_commondata_multi_index, std::string spat_pos, std::string mat_pos, bool newton_cotes=false)
 
int getRule (int order)
 
double area2D (double *a, double *b, double *c)
 
virtual MoFEMErrorCode setGaussPts (int order)
 
 ~MortarContactElement ()
 
- Public Member Functions inherited from SimpleContactProblem::ConvectMasterContactElement
 ConvectMasterContactElement (MoFEM::Interface &m_field, string spat_pos, string mat_pos, bool newton_cotes=false)
 
boost::shared_ptr< ConvectSlaveIntegrationPts > getConvectPtr ()
 
int getRule (int order)
 
MoFEMErrorCode setGaussPts (int order)
 
- Public Member Functions inherited from SimpleContactProblem::SimpleContactElement
 SimpleContactElement (MoFEM::Interface &m_field, bool newton_cotes=false)
 
MoFEMErrorCode preProcess ()
 
MoFEMErrorCode postProcess ()
 
int getRule (int order)
 

Public Attributes

MoFEM::Interface & mField
 
boost::shared_ptr< ContactSearchKdTree::ContactCommonData_multiIndex > contactCommondataMultiIndex
 
bool newtonCotes
 
- Public Attributes inherited from SimpleContactProblem::SimpleContactElement
MoFEM::Interface & mField
 
bool newtonCotes
 
SmartPetscObj< Vec > contactStateVec
 
friend ConvectSlaveIntegrationPts
 

Additional Inherited Members

- Protected Attributes inherited from SimpleContactProblem::ConvectMasterContactElement
boost::shared_ptr< ConvectSlaveIntegrationPts > convectPtr
 

Detailed Description

Definition at line 29 of file MortarContactProblem.hpp.

Constructor & Destructor Documentation

◆ MortarContactElement()

MortarContactProblem::MortarContactElement::MortarContactElement ( MoFEM::Interface &  m_field,
boost::shared_ptr< ContactSearchKdTree::ContactCommonData_multiIndex >  contact_commondata_multi_index,
std::string  spat_pos,
std::string  mat_pos,
bool  newton_cotes = false 
)
inline

Definition at line 36 of file MortarContactProblem.hpp.

41  : SimpleContactProblem::ConvectMasterContactElement(
42  m_field, spat_pos, mat_pos, newton_cotes),
43  mField(m_field),
44  contactCommondataMultiIndex(contact_commondata_multi_index),
45  newtonCotes(newton_cotes) {}

◆ ~MortarContactElement()

MortarContactProblem::MortarContactElement::~MortarContactElement ( )
inline

Definition at line 58 of file MortarContactProblem.hpp.

58 {}

Member Function Documentation

◆ area2D()

double MortarContactProblem::MortarContactElement::area2D ( double *  a,
double *  b,
double *  c 
)
inline

Definition at line 50 of file MortarContactProblem.hpp.

50  {
51  // (b-a)x(c-a) / 2
52  return ((b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])) /
53  2;
54  }

◆ getRule()

int MortarContactProblem::MortarContactElement::getRule ( int  order)
inline

Definition at line 47 of file MortarContactProblem.hpp.

47 { return -1; };

◆ setGaussPts()

MoFEMErrorCode MortarContactProblem::MortarContactElement::setGaussPts ( int  order)
virtual

Reimplemented from SimpleContactProblem::SimpleContactElement.

Reimplemented in MortarContactProblem::MortarConvectSlaveContactElement, and MortarContactProblem::MortarConvectMasterContactElement.

Definition at line 16 of file MortarContactProblem.cpp.

16  {
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 }

Member Data Documentation

◆ contactCommondataMultiIndex

boost::shared_ptr<ContactSearchKdTree::ContactCommonData_multiIndex> MortarContactProblem::MortarContactElement::contactCommondataMultiIndex

Definition at line 33 of file MortarContactProblem.hpp.

◆ mField

MoFEM::Interface& MortarContactProblem::MortarContactElement::mField

Definition at line 31 of file MortarContactProblem.hpp.

◆ newtonCotes

bool MortarContactProblem::MortarContactElement::newtonCotes

Definition at line 35 of file MortarContactProblem.hpp.

The documentation for this struct was generated from the following files:
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
CHKERRQ_MOAB
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:454
NOBASE
@ NOBASE
Definition: definitions.h:59
QUAD_2D_TABLE
static QUAD *const QUAD_2D_TABLE[]
Definition: quad.h:175
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
ContactSearchKdTree
Definition: ContactSearchKdTree.hpp:18
MortarContactProblem::MortarContactElement::area2D
double area2D(double *a, double *b, double *c)
Definition: MortarContactProblem.hpp:50
FTensor::Tensor2< double *, 3, 3 >
order
constexpr int order
Definition: dg_projection.cpp:18
ShapeDiffMBTRI
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:194
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MortarContactProblem::MortarContactElement::newtonCotes
bool newtonCotes
Definition: MortarContactProblem.hpp:35
a
constexpr double a
Definition: approx_sphere.cpp:30
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
ContactSearchKdTree::Prism_tag
Definition: ContactSearchKdTree.hpp:44
QUAD_::npoints
int npoints
Definition: quad.h:29
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
FTensor::Index< 'i', 3 >
MortarContactProblem::MortarContactElement::contactCommondataMultiIndex
boost::shared_ptr< ContactSearchKdTree::ContactCommonData_multiIndex > contactCommondataMultiIndex
Definition: MortarContactProblem.hpp:33
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
Range
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
SimpleContactProblem::ConvectMasterContactElement
Element used to integrate on slave surfaces. It convects integration points on slaves,...
Definition: SimpleContact.hpp:181
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
MortarContactProblem::MortarContactElement::mField
MoFEM::Interface & mField
Definition: MortarContactProblem.hpp:31
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
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
Generated by Doxygen 1.8.17 and hosted at University of Glasgow