v0.14.0
Loading...
Searching...
No Matches
ContactSearchKdTree.hpp
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#ifndef __CONTACT_SEARCH_KD_TREE__
16#define __CONTACT_SEARCH_KD_TREE__
17
19
22 : mField(m_field), kdTree(&m_field.get_moab()) {}
23
25 AdaptiveKDTree kdTree;
27 // kd-tree from
28 MoFEMErrorCode buildTree(Range &range_surf_master) {
30 CHKERR kdTree.build_tree(range_surf_master, &kdTreeRootMeshset);
32 }
33
34 // Definitition of multi-index containers consisting of prisms between master
35 // and slave trias and their corresponding triangles integration_tris
39 ContactCommonData(EntityHandle _prism, Range _common_integrated_triangle)
40 : pRism(_prism), commonIntegratedTriangle(_common_integrated_triangle) {
41 }
42 };
43
44 struct Prism_tag {};
45
46 typedef multi_index_container<
47 boost::shared_ptr<ContactCommonData>,
48 indexed_by<
49 ordered_unique<tag<Prism_tag>, member<ContactCommonData, EntityHandle,
52
53 // Rotation matrix to transform the 3D triangle to 2D triangle (because
54 // clipper library only works in 2D)
55 MoFEMErrorCode rotationMatrix(MatrixDouble &m_rot, VectorDouble &v_coords) {
57
61
62 VectorDouble3 v_s1(3), v_s2(3), v_n(3);
63 double diffN[6];
65
66 CHKERR ShapeFaceBaseMBTRI(diffN, &*v_coords.data().begin(),
67 &*v_n.data().begin(), &*v_s1.data().begin(),
68 &*v_s2.data().begin());
69
70 FTensor::Tensor1<double *, 3> t_s1(&v_s1[0], &v_s1[1], &v_s1[2], 3);
71 FTensor::Tensor1<double *, 3> t_s2(&v_s2[0], &v_s2[1], &v_s2[2], 3);
72 FTensor::Tensor1<double *, 3> t_n(&v_n[0], &v_n[1], &v_n[2], 3);
73
74 double s1_mag = sqrt(t_s1(i) * t_s1(i));
75 double n_mag = sqrt(t_n(i) * t_n(i));
76 v_s1 /= s1_mag;
77 v_n /= n_mag;
78
79 t_s2(i) = FTensor::levi_civita(i, j, k) * t_n(j) * t_s1(k);
80
81 for (int ii = 0; ii != 3; ++ii) {
82 m_rot(0, ii) = t_s1(ii);
83 m_rot(1, ii) = t_s2(ii);
84 m_rot(2, ii) = t_n(ii);
85 }
86
88 }
89
90 MoFEMErrorCode getNormalVector(VectorDouble &n_vec, VectorDouble &v_coords) {
92
94 double s1[3], s2[3];
95 double diffN[6];
97
98 CHKERR ShapeFaceBaseMBTRI(diffN, &*v_coords.data().begin(),
99 &*n_vec.data().begin(), s1, s2);
100
101 FTensor::Tensor1<double *, 3> t_n(&n_vec[0], &n_vec[1], &n_vec[2], 3);
102 double n_mag = sqrt(t_n(i) * t_n(i));
103 n_vec /= n_mag;
104
106 }
107
108 // Searching for contacting triangles with the help of kdtree from MOAB
109 // fill the multi-index container (contact_commondata_multi_index) conissting
110 // of prism with corresponding integration triangles fill
111 // range_slave_master_prisms consisting of all the prisms (to be used for
112 // looping over all these finite elements)
113 PetscErrorCode
114 contactSearchAlgorithm(Range &range_surf_master, Range &range_surf_slave,
115 boost::shared_ptr<ContactCommonData_multiIndex>
116 contact_commondata_multi_index,
117 Range &range_slave_master_prisms) {
119
122 auto get_tensor_vec = [](VectorDouble &v) {
124 &v(2));
125 };
126
127 // clipper library way to define one master and one slave tri and the out
128 // put will store in solution
129 Paths path_slave(1), path_master(1), solution;
130
131 // it is used to convert double to integer to be used in clipper library (as
132 // the library only works for integer) so we will convert all dimensions
133 // from double to integer and after finding common polygon will convert back
134 // to double
135 double roundfact = 1e6;
136
137 // Ranges of polygons in each intersections and created tris
138 Range tris_in_polygon;
139
140 VectorDouble v_slave_tri_cent;
141 v_slave_tri_cent.resize(3, false);
142
143 VectorDouble v_coords_slave_new;
144 v_coords_slave_new.resize(9, false);
145
146 MatrixDouble m_rot;
147 m_rot.resize(3, 3, false);
148
149 VectorDouble n_vec;
150 n_vec.resize(3, false);
151
152 VectorDouble v_coords_slave;
153 v_coords_slave.resize(9, false);
154
155 // loop over slave triagles (to create prisms)
156 int slave_num = 0;
157 for (Range::iterator tri_it_slave = range_surf_slave.begin();
158 tri_it_slave != range_surf_slave.end(); ++tri_it_slave) {
159 ++slave_num;
160
161 int num_nodes_slave;
162 const EntityHandle *conn_slave = NULL;
163 CHKERR mField.get_moab().get_connectivity(*tri_it_slave, conn_slave,
164 num_nodes_slave);
165
166 // insert prisms between master and slave triangles, who have common
167 // polygons last three nodes belong to the slave tri and the other three
168 // belong to the master tris
169 EntityHandle slave_master_1prism;
170 EntityHandle prism_nodes[6];
171 // slave is the top tri so add it to the top tri of prism [ii+3]
172 for (int ii = 0; ii != 3; ++ii) {
173 prism_nodes[ii + 3] = conn_slave[ii];
174 }
175
176 CHKERR mField.get_moab().get_coords(conn_slave, 3,
177 &*v_coords_slave.data().begin());
178
179 auto t_coords_slave = get_tensor_vec(v_coords_slave);
180 // rotation matrix (3x3) which will be used to transform the tri to z
181 // plane (this is the same for both slave and master tris) so
182 // calculate it only once
183 CHKERR rotationMatrix(m_rot, v_coords_slave);
185 &m_rot(0, 0), &m_rot(0, 1), &m_rot(0, 2), &m_rot(1, 0), &m_rot(1, 1),
186 &m_rot(1, 2), &m_rot(2, 0), &m_rot(2, 1), &m_rot(2, 2)};
187
188 CHKERR getNormalVector(n_vec, v_coords_slave);
189 auto t_slave_normal = get_tensor_vec(n_vec);
190
191 // as we need to transfer the coordinate to z=0
192 double z_shift;
193
194 auto t_coords_slave_new = get_tensor_vec(v_coords_slave_new);
195
196 auto t_slave_tri_cent = get_tensor_vec(v_slave_tri_cent);
197 t_slave_tri_cent(i) = 0;
198
199 // transfer the slave tri coordinates (n=3) one by one
200 for (int ii = 0; ii != 3; ++ii) {
201 t_coords_slave_new(i) = t_m_rot(i, j) * t_coords_slave(j);
202 t_slave_tri_cent(i) += t_coords_slave(i) / 3;
203 if (ii == 0)
204 z_shift = v_coords_slave_new[2];
205
206 t_coords_slave_new(2) = 0;
207 ++t_coords_slave_new;
208 ++t_coords_slave;
209 }
210
211 // multiply with roundfact to convert to integer
212 v_coords_slave_new = roundfact * v_coords_slave_new;
213
214 // clear clipper vector for slave tri and fill it to be used in clipper
215 path_slave[0].clear();
216 path_slave[0] << IntPoint(v_coords_slave_new[0], v_coords_slave_new[1])
217 << IntPoint(v_coords_slave_new[3], v_coords_slave_new[4])
218 << IntPoint(v_coords_slave_new[6], v_coords_slave_new[7]);
219
220 // find centroid of triangle based on original coordinates (this is
221 // created to search the nearest master tri) it is to be used in kd-tree
222
223 // this range will be filled in recursive manner, starting from the
224 // nearset one and unitill there is no more intersecting master tris
225 Range range_closest_master_tris;
226 range_closest_master_tris.clear();
227 double closest_point[3];
228 EntityHandle closest_master_tri;
229 CHKERR kdTree.closest_triangle(kdTreeRootMeshset,
230 &*v_slave_tri_cent.data().begin(),
231 closest_point, closest_master_tri);
232
233 range_closest_master_tris.insert(closest_master_tri);
234
235 // this while loop will continue until there is no more master tri for
236 // slave tris
237 bool flag_first_master =
238 true; // for the first master, we will not find the adjacenency
239 int master_num = 0;
240 bool flag_end_of_search = false;
241 bool flag_temp = true;
242 while (flag_end_of_search == false) {
243
244 Range range_master_tris_on_surf_new;
245 if (flag_first_master == true) { // first master
246 range_master_tris_on_surf_new = range_closest_master_tris;
247 flag_first_master = false;
248 } else { // the rest of master determine from adjacencies
249 flag_temp = false;
250 Range range_conn;
251 CHKERR mField.get_moab().get_connectivity(range_closest_master_tris,
252 range_conn);
253
254 Range range_adj_entities;
255 // in the get_adjacencies, if moab::Interface::INTERSECT it will give
256 // only the tri, which is interesect of the adjacency of all nodes
257 // moab::Interface::UNION will give all the tri, even the one which
258 // are not on the surface, so then need to intersect with the master
259 // to get what we want
260 CHKERR mField.get_moab().get_adjacencies(
261 range_conn, 2, false, range_adj_entities, moab::Interface::UNION);
262
263 // triangle on the surface, i.e. removing the tri inside the volume
264 Range range_master_tris_on_surf =
265 intersect(range_surf_master, range_adj_entities);
266
267 // removing the tri from the range_tri_on_surf, which already exists
268 // in range_closest_master_tris
269 range_master_tris_on_surf_new =
270 subtract(range_master_tris_on_surf, range_closest_master_tris);
271 }
272
273 // all used in master loop
274 VectorDouble v_coords_master;
275 v_coords_master.resize(9, false);
276
277 VectorDouble n_vec_master;
278 n_vec_master.resize(3, false);
279
280 // fill the master clipper vector
281 VectorDouble v_coords_master_new;
282 v_coords_master_new.resize(9, false);
283
284 double coord[3];
285 // looping over all master range_master_tris_on_surf_new and discart
286 // those elements, which are not in contact with master (will do it with
287 // clipper library)
288 int num_intersections = 0;
289 for (Range::iterator tri_it_master =
290 range_master_tris_on_surf_new.begin();
291 tri_it_master != range_master_tris_on_surf_new.end();
292 ++tri_it_master) {
293 const EntityHandle *conn_master = NULL;
294 int num_nodes_master = 0;
295 CHKERR mField.get_moab().get_connectivity(*tri_it_master, conn_master,
296 num_nodes_master);
297
298 CHKERR mField.get_moab().get_coords(conn_master, 3,
299 &*v_coords_master.data().begin());
300
301 auto t_coords_master = get_tensor_vec(v_coords_master);
302
303 CHKERR getNormalVector(n_vec_master, v_coords_master);
304
305 auto t_master_normal = get_tensor_vec(n_vec_master);
306
307 // 0.035 angle of 88 deg
308 if (fabs(t_slave_normal(i) * t_master_normal(i)) < 0.035)
309 continue;
310
311 auto t_coords_master_new = get_tensor_vec(v_coords_master_new);
312
313 // rotation matrix for this master element is the same as for the
314 // above slave element
315 for (int ii = 0; ii != 3; ++ii) {
316 t_coords_master_new(i) = t_m_rot(i, j) * t_coords_master(j);
317 t_coords_master_new(2) = 0;
318 ++t_coords_master_new;
319 ++t_coords_master;
320 }
321
322 v_coords_master_new = roundfact * v_coords_master_new;
323 path_master[0].clear();
324 path_master[0]
325 << IntPoint(v_coords_master_new[0], v_coords_master_new[1])
326 << IntPoint(v_coords_master_new[3], v_coords_master_new[4])
327 << IntPoint(v_coords_master_new[6], v_coords_master_new[7]);
328
329 // perform intersection in clipper
330 Clipper c;
331 c.AddPaths(path_master, ptSubject, true);
332 c.AddPaths(path_slave, ptClip, true);
333 c.Execute(ctIntersection, solution, pftNonZero, pftNonZero);
334
335 if (solution.size() != 1 && flag_temp) {
336 double ray_dir[3];
337 CHKERR Tools::getTriNormal(&v_coords_slave(0), ray_dir);
338 double norm = cblas_dnrm2(3, ray_dir, 1);
339 cblas_dscal(3, 1. / norm, ray_dir, 1);
340 const double tol = 1e-6;
341 std::vector<double> distance;
342 std::vector<EntityHandle> triangles_out;
343 CHKERR kdTree.ray_intersect_triangles(kdTreeRootMeshset, tol, ray_dir,
344 &*v_slave_tri_cent.data().begin(),
345 triangles_out, distance);
346
347 if (!triangles_out.empty()) {
348 range_closest_master_tris.clear();
349 range_closest_master_tris.insert(triangles_out.front());
350 flag_first_master = true;
351 num_intersections = -1;
352 break;
353 }
354 }
355
356 if (solution.size() == 1) { // if it intersect
357 range_closest_master_tris.insert(
358 *tri_it_master); // update range_closest_master_tris
359 num_intersections += 1;
360
361 int n_vertices_polygon = solution[0].size();
362 // the first three nodes for the prism are the master tri (insert
363 // prisms only if there is intersection between slave and master
364 // tris)
365 for (int ii = 0; ii != 3; ++ii) {
366 prism_nodes[ii] = conn_master[ii];
367 }
368
369 // creating prism between master and slave tris
370 CHKERR mField.get_moab().create_element(MBPRISM, prism_nodes, 6,
371 slave_master_1prism);
372
373 range_slave_master_prisms.insert(slave_master_1prism);
374
375 // coordinates of intersection polygon between master and slave
376 // here we have only coordinates of polygon, so we will need to
377 // construct vertices and polygon from it
378 VectorDouble v_coord_polygon;
379 v_coord_polygon.resize(3 * n_vertices_polygon);
380
381 // perform transfermation (tansfer all coordinates of polygon)
382 // create moab vertices with given coordinates (from clippper) and
383 // then create polygon from these vertices
384 VectorDouble v_coord_polygon_new;
385 v_coord_polygon_new.resize(3 * n_vertices_polygon);
386
387 auto t_coord_polygon_new = get_tensor_vec(v_coord_polygon_new);
388 auto t_coord_polygon = get_tensor_vec(v_coord_polygon);
389
390 std::vector<EntityHandle> conn_polygon(n_vertices_polygon);
391 for (int ii = 0; ii != n_vertices_polygon; ++ii) {
392 t_coord_polygon(0) = solution[0][ii].X / roundfact;
393 t_coord_polygon(1) = solution[0][ii].Y / roundfact;
394 t_coord_polygon(2) += z_shift;
395
396 t_coord_polygon_new(i) = t_m_rot(j, i) * t_coord_polygon(j);
397
398 coord[0] = t_coord_polygon_new(0);
399 coord[1] = t_coord_polygon_new(1);
400 coord[2] = t_coord_polygon_new(2);
401 CHKERR mField.get_moab().create_vertex(coord, conn_polygon[ii]);
402
403 ++t_coord_polygon_new;
404 ++t_coord_polygon;
405 }
406
407 EntityHandle polygon;
408 CHKERR mField.get_moab().create_element(
409 MBPOLYGON, &conn_polygon[0], n_vertices_polygon, polygon);
410
411 // triangulate the polygon (if it only have 3 vertices, so it is
412 // triangle, so don't need to triangulate)
413 tris_in_polygon.clear();
414 if (n_vertices_polygon == 3) {
415 // tris_in_polygon.insert(polygon);
416 EntityHandle new_tri_conn[3];
417 new_tri_conn[0] = conn_polygon[0];
418 new_tri_conn[1] = conn_polygon[1];
419 new_tri_conn[2] = conn_polygon[2];
420 EntityHandle new_tri;
421 CHKERR mField.get_moab().create_element(MBTRI, new_tri_conn, 3,
422 new_tri);
423
424 tris_in_polygon.insert(new_tri);
425 } else {
426 EntityHandle new_tri_conn[3];
427 new_tri_conn[0] = conn_polygon[0];
428 for (int ii = 0; ii != n_vertices_polygon - 2; ++ii) {
429 new_tri_conn[1] = conn_polygon[ii + 1];
430 new_tri_conn[2] = conn_polygon[ii + 2];
431 // create tri in polygon
432 EntityHandle new_tri;
433 CHKERR mField.get_moab().create_element(MBTRI, new_tri_conn, 3,
434 new_tri);
435
436 tris_in_polygon.insert(new_tri);
437 }
438 }
439
440 // fill the multi-index container with pair (prism between slave and
441 // master, triangles in the common polygon)
442 contact_commondata_multi_index->insert(
443 boost::shared_ptr<ContactCommonData>(new ContactCommonData(
444 slave_master_1prism, tris_in_polygon)));
445 }
446 } // master loop for new layer arrond the previous layer
447
448 if (num_intersections == 0)
449 flag_end_of_search = true; // this will stop further searching for the
450 // current slave and will go to next one
451 } // while loop
452 } // slave loop
453
455 }
456};
457#endif
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
PetscErrorCode ShapeFaceBaseMBTRI(double *diffN, const double *coords, double *normal, double *s1, double *s2)
Definition: fem_tools.c:204
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:194
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
double tol
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
ContactCommonData(EntityHandle _prism, Range _common_integrated_triangle)
ContactSearchKdTree(MoFEM::Interface &m_field)
MoFEM::Interface & mField
PetscErrorCode contactSearchAlgorithm(Range &range_surf_master, Range &range_surf_slave, boost::shared_ptr< ContactCommonData_multiIndex > contact_commondata_multi_index, Range &range_slave_master_prisms)
MoFEMErrorCode buildTree(Range &range_surf_master)
MoFEMErrorCode rotationMatrix(MatrixDouble &m_rot, VectorDouble &v_coords)
multi_index_container< boost::shared_ptr< ContactCommonData >, indexed_by< ordered_unique< tag< Prism_tag >, member< ContactCommonData, EntityHandle, &ContactCommonData::pRism > > > > ContactCommonData_multiIndex
MoFEMErrorCode getNormalVector(VectorDouble &n_vec, VectorDouble &v_coords)
virtual moab::Interface & get_moab()=0
Deprecated interface functions.