Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "XFEMCutElem2D.h"
11 :
12 : #include "EFANode.h"
13 : #include "EFAEdge.h"
14 : #include "EFAFragment2D.h"
15 : #include "XFEMFuncs.h"
16 : #include "MooseError.h"
17 :
18 : #include "libmesh/mesh.h"
19 : #include "libmesh/elem.h"
20 : #include "libmesh/node.h"
21 : #include "libmesh/petsc_macro.h"
22 : #include "petscblaslapack.h"
23 :
24 9298 : XFEMCutElem2D::XFEMCutElem2D(Elem * elem,
25 : const EFAElement2D * const CEMelem,
26 : unsigned int n_qpoints,
27 9298 : unsigned int n_sides)
28 9298 : : XFEMCutElem(elem, n_qpoints, n_sides), _efa_elem2d(CEMelem, true)
29 : {
30 9298 : computePhysicalVolumeFraction();
31 9298 : }
32 :
33 18596 : XFEMCutElem2D::~XFEMCutElem2D() {}
34 :
35 : Point
36 4227238 : XFEMCutElem2D::getNodeCoordinates(EFANode * CEMnode, MeshBase * displaced_mesh) const
37 : {
38 : Point node_coor(0.0, 0.0, 0.0);
39 : std::vector<EFANode *> master_nodes;
40 : std::vector<Point> master_points;
41 : std::vector<double> master_weights;
42 :
43 4227238 : _efa_elem2d.getMasterInfo(CEMnode, master_nodes, master_weights);
44 12321762 : for (unsigned int i = 0; i < master_nodes.size(); ++i)
45 : {
46 8094524 : if (master_nodes[i]->category() == EFANode::N_CATEGORY_LOCAL_INDEX)
47 : {
48 8094524 : Node * node = _nodes[master_nodes[i]->id()];
49 8094524 : if (displaced_mesh)
50 5719292 : node = displaced_mesh->node_ptr(node->id());
51 8094524 : Point node_p((*node)(0), (*node)(1), 0.0);
52 8094524 : master_points.push_back(node_p);
53 : }
54 : else
55 0 : mooseError("master nodes must be local");
56 : }
57 12321762 : for (unsigned int i = 0; i < master_nodes.size(); ++i)
58 : node_coor += master_weights[i] * master_points[i];
59 4227238 : return node_coor;
60 4227238 : }
61 :
62 : void
63 92665 : XFEMCutElem2D::computePhysicalVolumeFraction()
64 : {
65 : Real frag_vol = 0.0;
66 :
67 : // Calculate area of entire element and fragment using the formula:
68 : // A = 1/2 sum_{i=0}^{n-1} (x_i y_{i+1} - x_{i+1} y{i})
69 :
70 453274 : for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
71 : {
72 360609 : Point edge_p1 = getNodeCoordinates(_efa_elem2d.getFragmentEdge(0, i)->getNode(0));
73 360609 : Point edge_p2 = getNodeCoordinates(_efa_elem2d.getFragmentEdge(0, i)->getNode(1));
74 360609 : frag_vol += 0.5 * (edge_p1(0) - edge_p2(0)) * (edge_p1(1) + edge_p2(1));
75 : }
76 92665 : _physical_volfrac = frag_vol / _elem_volume;
77 92665 : }
78 :
79 : void
80 18 : XFEMCutElem2D::computePhysicalFaceAreaFraction(unsigned int side)
81 : {
82 : Real frag_surf = 0.0;
83 :
84 18 : EFAEdge * edge = _efa_elem2d.getEdge(side);
85 :
86 42 : for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
87 : {
88 42 : EFANode * node_1 = _efa_elem2d.getFragmentEdge(0, i)->getNode(0);
89 42 : EFANode * node_2 = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
90 :
91 : /// find a fragment edge which is covered by element side
92 42 : if (edge->containsNode(node_1) && edge->containsNode(node_2))
93 : {
94 18 : Point edge_p1 = getNodeCoordinates(node_1);
95 18 : Point edge_p2 = getNodeCoordinates(node_2);
96 18 : frag_surf = (edge_p1 - edge_p2).norm();
97 18 : _physical_areafrac[side] = frag_surf / _elem_side_area[side];
98 : return;
99 : }
100 : }
101 0 : _physical_areafrac[side] = 1.0;
102 : }
103 :
104 : void
105 348 : XFEMCutElem2D::computeMomentFittingWeights()
106 : {
107 : // Purpose: calculate new weights via moment-fitting method
108 348 : std::vector<Point> elem_nodes(_n_nodes, Point(0.0, 0.0, 0.0));
109 : std::vector<std::vector<Real>> wsg;
110 :
111 1740 : for (unsigned int i = 0; i < _n_nodes; ++i)
112 1392 : elem_nodes[i] = (*_nodes[i]);
113 :
114 348 : if (_efa_elem2d.isPartial() && _n_qpoints <= 6) // ONLY work for <= 6 q_points
115 : {
116 : std::vector<std::vector<Real>> tsg;
117 348 : getPhysicalQuadraturePoints(tsg); // get tsg - QPs within partial element
118 348 : solveMomentFitting(
119 : _n_nodes, _n_qpoints, elem_nodes, tsg, wsg); // get wsg - QPs from moment-fitting
120 348 : _new_weights.resize(wsg.size(), 1.0);
121 1776 : for (unsigned int i = 0; i < wsg.size(); ++i)
122 1428 : _new_weights[i] = wsg[i][2]; // weight multiplier
123 348 : }
124 : else
125 0 : _new_weights.resize(_n_qpoints, _physical_volfrac);
126 348 : }
127 :
128 : Point
129 486576 : XFEMCutElem2D::getCutPlaneOrigin(unsigned int plane_id, MeshBase * displaced_mesh) const
130 : {
131 : Point orig(0.0, 0.0, 0.0);
132 : std::vector<std::vector<EFANode *>> cut_line_nodes;
133 2387062 : for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
134 : {
135 1900486 : if (_efa_elem2d.getFragment(0)->isEdgeInterior(i))
136 : {
137 490644 : std::vector<EFANode *> node_line(2, nullptr);
138 490644 : node_line[0] = _efa_elem2d.getFragmentEdge(0, i)->getNode(0);
139 490644 : node_line[1] = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
140 490644 : cut_line_nodes.push_back(node_line);
141 490644 : }
142 : }
143 486576 : if (cut_line_nodes.size() == 0)
144 0 : mooseError("no cut line found in this element");
145 486576 : if (plane_id < cut_line_nodes.size()) // valid plane_id
146 331254 : orig = getNodeCoordinates(cut_line_nodes[plane_id][0], displaced_mesh);
147 486576 : return orig;
148 486576 : }
149 :
150 : Point
151 1113278 : XFEMCutElem2D::getCutPlaneNormal(unsigned int plane_id, MeshBase * displaced_mesh) const
152 : {
153 : Point normal(0.0, 0.0, 0.0);
154 : std::vector<std::vector<EFANode *>> cut_line_nodes;
155 5500892 : for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
156 : {
157 4387614 : if (_efa_elem2d.getFragment(0)->isEdgeInterior(i))
158 : {
159 1117346 : std::vector<EFANode *> node_line(2, nullptr);
160 1117346 : node_line[0] = _efa_elem2d.getFragmentEdge(0, i)->getNode(0);
161 1117346 : node_line[1] = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
162 1117346 : cut_line_nodes.push_back(node_line);
163 1117346 : }
164 : }
165 1113278 : if (cut_line_nodes.size() == 0)
166 0 : mooseError("no cut line found in this element");
167 1113278 : if (plane_id < cut_line_nodes.size()) // valid plane_id
168 : {
169 957956 : Point cut_line_p1 = getNodeCoordinates(cut_line_nodes[plane_id][0], displaced_mesh);
170 957956 : Point cut_line_p2 = getNodeCoordinates(cut_line_nodes[plane_id][1], displaced_mesh);
171 : Point cut_line = cut_line_p2 - cut_line_p1;
172 957956 : Real len = std::sqrt(cut_line.norm_sq());
173 : cut_line /= len;
174 957956 : normal = Point(cut_line(1), -cut_line(0), 0.0);
175 : }
176 1113278 : return normal;
177 1113278 : }
178 :
179 : void
180 2011 : XFEMCutElem2D::getCrackTipOriginAndDirection(unsigned tip_id,
181 : Point & origin,
182 : Point & direction) const
183 : {
184 : // TODO: two cut plane case is not working
185 : std::vector<EFANode *> cut_line_nodes;
186 10033 : for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
187 : {
188 8022 : if (_efa_elem2d.getFragment(0)->isEdgeInterior(i))
189 : {
190 2011 : std::vector<EFANode *> node_line(2, nullptr);
191 2011 : node_line[0] = _efa_elem2d.getFragmentEdge(0, i)->getNode(0);
192 2011 : node_line[1] = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
193 2011 : if (node_line[1]->id() == tip_id)
194 : {
195 987 : cut_line_nodes.push_back(node_line[0]);
196 987 : cut_line_nodes.push_back(node_line[1]);
197 : }
198 1024 : else if (node_line[0]->id() == tip_id)
199 : {
200 1024 : node_line[1] = node_line[0];
201 1024 : node_line[0] = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
202 1024 : cut_line_nodes.push_back(node_line[0]);
203 1024 : cut_line_nodes.push_back(node_line[1]);
204 : }
205 2011 : }
206 : }
207 2011 : if (cut_line_nodes.size() == 0)
208 0 : mooseError("no cut line found in this element");
209 :
210 2011 : Point cut_line_p1 = getNodeCoordinates(cut_line_nodes[0]);
211 2011 : Point cut_line_p2 = getNodeCoordinates(cut_line_nodes[1]);
212 : Point cut_line = cut_line_p2 - cut_line_p1;
213 2011 : Real len = std::sqrt(cut_line.norm_sq());
214 : cut_line /= len;
215 2011 : origin = cut_line_p2;
216 2011 : direction = Point(cut_line(0), cut_line(1), 0.0);
217 2011 : }
218 :
219 : void
220 0 : XFEMCutElem2D::getFragmentFaces(std::vector<std::vector<Point>> & frag_faces,
221 : MeshBase * displaced_mesh) const
222 : {
223 0 : frag_faces.clear();
224 0 : for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
225 : {
226 0 : std::vector<Point> edge_points(2, Point(0.0, 0.0, 0.0));
227 0 : edge_points[0] =
228 0 : getNodeCoordinates(_efa_elem2d.getFragmentEdge(0, i)->getNode(0), displaced_mesh);
229 0 : edge_points[1] =
230 0 : getNodeCoordinates(_efa_elem2d.getFragmentEdge(0, i)->getNode(1), displaced_mesh);
231 0 : frag_faces.push_back(edge_points);
232 0 : }
233 0 : }
234 :
235 : const EFAElement *
236 1298225 : XFEMCutElem2D::getEFAElement() const
237 : {
238 1298225 : return &_efa_elem2d;
239 : }
240 :
241 : unsigned int
242 171864 : XFEMCutElem2D::numCutPlanes() const
243 : {
244 : unsigned int counter = 0;
245 859426 : for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
246 687562 : if (_efa_elem2d.getFragment(0)->isEdgeInterior(i))
247 171864 : counter += 1;
248 171864 : return counter;
249 : }
250 :
251 : void
252 348 : XFEMCutElem2D::getPhysicalQuadraturePoints(std::vector<std::vector<Real>> & tsg)
253 : {
254 : // Get the coords for parial element nodes
255 348 : EFAFragment2D * frag = _efa_elem2d.getFragment(0);
256 348 : unsigned int nnd_pe = frag->numEdges();
257 348 : std::vector<Point> frag_points(nnd_pe, Point(0.0, 0.0, 0.0)); // nodal coord of partial elem
258 348 : Real jac = 0.0;
259 :
260 1740 : for (unsigned int j = 0; j < nnd_pe; ++j)
261 1392 : frag_points[j] = getNodeCoordinates(frag->getEdge(j)->getNode(0));
262 :
263 : // Get centroid coords for partial elements
264 : Point xcrd(0.0, 0.0, 0.0);
265 1740 : for (unsigned int j = 0; j < nnd_pe; ++j)
266 1392 : xcrd += frag_points[j];
267 348 : xcrd *= (1.0 / nnd_pe);
268 :
269 : // Get tsg - the physical coords of Gaussian Q-points for partial elements
270 348 : if ((nnd_pe == 3) || (nnd_pe == 4)) // partial element is a triangle or quad
271 : {
272 : std::vector<std::vector<Real>> sg2;
273 276 : std::vector<std::vector<Real>> shape(nnd_pe, std::vector<Real>(3, 0.0));
274 276 : Xfem::stdQuadr2D(nnd_pe, 2, sg2);
275 1308 : for (unsigned int l = 0; l < sg2.size(); ++l)
276 : {
277 1032 : Xfem::shapeFunc2D(nnd_pe, sg2[l], frag_points, shape, jac, true); // Get shape
278 1032 : std::vector<Real> tsg_line(3, 0.0);
279 4944 : for (unsigned int k = 0; k < nnd_pe; ++k)
280 : {
281 3912 : tsg_line[0] += shape[k][2] * frag_points[k](0);
282 3912 : tsg_line[1] += shape[k][2] * frag_points[k](1);
283 : }
284 1032 : if (nnd_pe == 3) // tri partial elem
285 216 : tsg_line[2] = sg2[l][3] * jac; // total weights
286 : else // quad partial elem
287 816 : tsg_line[2] = sg2[l][2] * jac; // total weights
288 1032 : tsg.push_back(tsg_line);
289 1032 : }
290 276 : }
291 72 : else if (nnd_pe >= 5) // partial element is a polygon
292 : {
293 432 : for (unsigned int j = 0; j < nnd_pe; ++j) // loop all sub-tris
294 : {
295 360 : std::vector<std::vector<Real>> shape(3, std::vector<Real>(3, 0.0));
296 360 : std::vector<Point> subtri_points(3, Point(0.0, 0.0, 0.0)); // sub-tri nodal coords
297 :
298 360 : int jplus1 = j < nnd_pe - 1 ? j + 1 : 0;
299 360 : subtri_points[0] = xcrd;
300 360 : subtri_points[1] = frag_points[j];
301 360 : subtri_points[2] = frag_points[jplus1];
302 :
303 : std::vector<std::vector<Real>> sg2;
304 360 : Xfem::stdQuadr2D(3, 2, sg2); // get sg2
305 1440 : for (unsigned int l = 0; l < sg2.size(); ++l) // loop all int pts on a sub-tri
306 : {
307 1080 : Xfem::shapeFunc2D(3, sg2[l], subtri_points, shape, jac, true); // Get shape
308 1080 : std::vector<Real> tsg_line(3, 0.0);
309 4320 : for (unsigned int k = 0; k < 3; ++k) // loop sub-tri nodes
310 : {
311 3240 : tsg_line[0] += shape[k][2] * subtri_points[k](0);
312 3240 : tsg_line[1] += shape[k][2] * subtri_points[k](1);
313 : }
314 1080 : tsg_line[2] = sg2[l][3] * jac;
315 1080 : tsg.push_back(tsg_line);
316 1080 : }
317 360 : }
318 : }
319 : else
320 0 : mooseError("Invalid partial element!");
321 348 : }
322 :
323 : void
324 348 : XFEMCutElem2D::solveMomentFitting(unsigned int nen,
325 : unsigned int nqp,
326 : std::vector<Point> & elem_nodes,
327 : std::vector<std::vector<Real>> & tsg,
328 : std::vector<std::vector<Real>> & wsg)
329 : {
330 : // Get physical coords for the new six-point rule
331 348 : std::vector<std::vector<Real>> shape(nen, std::vector<Real>(3, 0.0));
332 : std::vector<std::vector<Real>> wss;
333 :
334 348 : if (nen == 4)
335 : {
336 348 : wss.resize(_qp_points.size());
337 1776 : for (unsigned int i = 0; i < _qp_points.size(); i++)
338 : {
339 1428 : wss[i].resize(3);
340 1428 : wss[i][0] = _qp_points[i](0);
341 1428 : wss[i][1] = _qp_points[i](1);
342 1428 : wss[i][2] = _qp_weights[i];
343 : }
344 : }
345 0 : else if (nen == 3)
346 : {
347 0 : wss.resize(_qp_points.size());
348 0 : for (unsigned int i = 0; i < _qp_points.size(); i++)
349 : {
350 0 : wss[i].resize(4);
351 0 : wss[i][0] = _qp_points[i](0);
352 0 : wss[i][1] = _qp_points[i](1);
353 0 : wss[i][2] = 1.0 - _qp_points[i](0) - _qp_points[i](1);
354 0 : wss[i][3] = _qp_weights[i];
355 : }
356 : }
357 : else
358 0 : mooseError("Invalid element");
359 :
360 348 : wsg.resize(wss.size());
361 1776 : for (unsigned int i = 0; i < wsg.size(); ++i)
362 1428 : wsg[i].resize(3, 0.0);
363 348 : Real jac = 0.0;
364 348 : std::vector<Real> old_weights(wss.size(), 0.0);
365 :
366 1776 : for (unsigned int l = 0; l < wsg.size(); ++l)
367 : {
368 1428 : Xfem::shapeFunc2D(nen, wss[l], elem_nodes, shape, jac, true); // Get shape
369 1428 : if (nen == 4) // 2D quad elem
370 1428 : old_weights[l] = wss[l][2] * jac; // weights for total element
371 : else if (nen == 3) // 2D triangle elem
372 0 : old_weights[l] = wss[l][3] * jac;
373 : else
374 : mooseError("Invalid element!");
375 7140 : for (unsigned int k = 0; k < nen; ++k) // physical coords of Q-pts
376 : {
377 5712 : wsg[l][0] += shape[k][2] * elem_nodes[k](0);
378 5712 : wsg[l][1] += shape[k][2] * elem_nodes[k](1);
379 : }
380 : }
381 :
382 : // Compute weights via moment fitting
383 : Real * A;
384 348 : A = new Real[wsg.size() * wsg.size()];
385 : unsigned ind = 0;
386 1776 : for (unsigned int i = 0; i < wsg.size(); ++i)
387 : {
388 1428 : A[ind] = 1.0; // const
389 1428 : if (nqp > 1)
390 1428 : A[1 + ind] = wsg[i][0]; // x
391 1428 : if (nqp > 2)
392 1428 : A[2 + ind] = wsg[i][1]; // y
393 1428 : if (nqp > 3)
394 1428 : A[3 + ind] = wsg[i][0] * wsg[i][1]; // x*y
395 1428 : if (nqp > 4)
396 108 : A[4 + ind] = wsg[i][0] * wsg[i][0]; // x^2
397 108 : if (nqp > 5)
398 108 : A[5 + ind] = wsg[i][1] * wsg[i][1]; // y^2
399 108 : if (nqp > 6)
400 0 : mooseError("Q-points of more than 6 are not allowed now!");
401 1428 : ind = ind + nqp;
402 : }
403 :
404 : Real * b;
405 348 : b = new Real[wsg.size()];
406 1776 : for (unsigned int i = 0; i < wsg.size(); ++i)
407 1428 : b[i] = 0.0;
408 2460 : for (unsigned int i = 0; i < tsg.size(); ++i)
409 : {
410 2112 : b[0] += tsg[i][2];
411 2112 : if (nqp > 1)
412 2112 : b[1] += tsg[i][2] * tsg[i][0];
413 2112 : if (nqp > 2)
414 2112 : b[2] += tsg[i][2] * tsg[i][1];
415 2112 : if (nqp > 3)
416 2112 : b[3] += tsg[i][2] * tsg[i][0] * tsg[i][1];
417 2112 : if (nqp > 4)
418 72 : b[4] += tsg[i][2] * tsg[i][0] * tsg[i][0];
419 72 : if (nqp > 5)
420 72 : b[5] += tsg[i][2] * tsg[i][1] * tsg[i][1];
421 72 : if (nqp > 6)
422 0 : mooseError("Q-points of more than 6 are not allowed now!");
423 : }
424 :
425 348 : PetscBLASInt nrhs = 1;
426 : PetscBLASInt info;
427 348 : PetscBLASInt n = wsg.size();
428 348 : std::vector<PetscBLASInt> ipiv(n);
429 :
430 348 : LAPACKgesv_(&n, &nrhs, A, &n, &ipiv[0], b, &n, &info);
431 :
432 1776 : for (unsigned int i = 0; i < wsg.size(); ++i)
433 1428 : wsg[i][2] = b[i] / old_weights[i]; // get the multiplier
434 :
435 : // delete arrays
436 348 : delete[] A;
437 348 : delete[] b;
438 348 : }
439 :
440 : void
441 626702 : XFEMCutElem2D::getIntersectionInfo(unsigned int plane_id,
442 : Point & normal,
443 : std::vector<Point> & intersectionPoints,
444 : MeshBase * displaced_mesh) const
445 : {
446 626702 : intersectionPoints.resize(2); // 2D: the intersection line is straight and can be represented by
447 : // two ending points.(may have issues with double cuts case!)
448 :
449 : std::vector<std::vector<EFANode *>> cut_line_nodes;
450 3113830 : for (unsigned int i = 0; i < _efa_elem2d.getFragment(0)->numEdges(); ++i)
451 : {
452 2487128 : if (_efa_elem2d.getFragment(0)->isEdgeInterior(i))
453 : {
454 626702 : std::vector<EFANode *> node_line(2, nullptr);
455 626702 : node_line[0] = _efa_elem2d.getFragmentEdge(0, i)->getNode(0);
456 626702 : node_line[1] = _efa_elem2d.getFragmentEdge(0, i)->getNode(1);
457 626702 : cut_line_nodes.push_back(node_line);
458 626702 : }
459 : }
460 626702 : if (cut_line_nodes.size() == 0)
461 0 : mooseError("No cut line found in this element");
462 :
463 626702 : if (plane_id < cut_line_nodes.size()) // valid plane_id
464 : {
465 626702 : intersectionPoints[0] = getNodeCoordinates(cut_line_nodes[plane_id][0], displaced_mesh);
466 626702 : intersectionPoints[1] = getNodeCoordinates(cut_line_nodes[plane_id][1], displaced_mesh);
467 : }
468 :
469 626702 : normal = getCutPlaneNormal(plane_id, displaced_mesh);
470 626702 : }
|
