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