www.mooseframework.org
XFEMCutElem.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 "XFEMCutElem.h"
9 
10 #include "libmesh/mesh.h"
11 #include "libmesh/elem.h"
12 #include "libmesh/node.h"
13 #include "libmesh/quadrature.h"
14 #include "libmesh/fe.h"
15 #include "libmesh/enum_fe_family.h"
16 
17 #include "EFANode.h"
18 #include "EFAElement.h"
19 
20 XFEMCutElem::XFEMCutElem(Elem * elem, unsigned int n_qpoints)
21  : _n_nodes(elem->n_nodes()), _n_qpoints(n_qpoints), _nodes(_n_nodes, NULL), _have_weights(false)
22 {
23  for (unsigned int i = 0; i < _n_nodes; ++i)
24  _nodes[i] = elem->get_node(i);
25  _elem_volume = elem->volume();
26 }
27 
29 
30 Real
32 {
33  return _physical_volfrac;
34 }
35 
36 void
37 XFEMCutElem::getWeightMultipliers(MooseArray<Real> & weights,
38  QBase * qrule,
39  Xfem::XFEM_QRULE xfem_qrule,
40  const MooseArray<Point> & q_points)
41 {
42  if (!_have_weights)
43  computeXFEMWeights(qrule, xfem_qrule, q_points);
44 
45  weights.resize(_new_weights.size());
46  for (unsigned int qp = 0; qp < _new_weights.size(); ++qp)
47  weights[qp] = _new_weights[qp];
48 }
49 
50 void
52  Xfem::XFEM_QRULE xfem_qrule,
53  const MooseArray<Point> & q_points)
54 {
55  _new_weights.clear();
56 
57  _new_weights.resize(qrule->n_points(), 1.0);
58 
59  switch (xfem_qrule)
60  {
61  case Xfem::VOLFRAC:
62  {
64  Real volfrac = getPhysicalVolumeFraction();
65  for (unsigned qp = 0; qp < qrule->n_points(); ++qp)
66  _new_weights[qp] = volfrac;
67  break;
68  }
70  {
71  // These are the coordinates in parametric coordinates
72  _qp_points = qrule->get_points();
73  _qp_weights = qrule->get_weights();
74 
76 
77  // Blend weights from moment fitting and volume fraction to avoid negative weights
78  Real alpha = 1.0;
79  if (*std::min_element(_new_weights.begin(), _new_weights.end()) < 0.0)
80  {
81  // One or more of the weights computed by moment fitting is negative.
82  // Blend moment and volume fraction weights to keep them nonnegative.
83  // Find the largest value of alpha that will keep all of the weights nonnegative.
84  for (unsigned int i = 0; i < _n_qpoints; ++i)
85  {
86  const Real denominator = _physical_volfrac - _new_weights[i];
87  if (denominator > 0.0) // Negative values would give a negative value for alpha, which
88  // must be between 0 and 1.
89  {
90  const Real alpha_i = _physical_volfrac / denominator;
91  if (alpha_i < alpha)
92  alpha = alpha_i;
93  }
94  }
95  for (unsigned int i = 0; i < _n_qpoints; ++i)
96  {
97  _new_weights[i] = alpha * _new_weights[i] + (1.0 - alpha) * _physical_volfrac;
98  if (_new_weights[i] < 0.0) // We can end up with small (roundoff) negative weights
99  _new_weights[i] = 0.0;
100  }
101  }
102  break;
103  }
104  case Xfem::DIRECT: // remove q-points outside the partial element's physical domain
105  {
106  bool nonzero = false;
107  for (unsigned qp = 0; qp < qrule->n_points(); ++qp)
108  {
109  // q_points contains quadrature point locations in physical coordinates
110  if (isPointPhysical(q_points[qp]))
111  {
112  nonzero = true;
113  _new_weights[qp] = 1.0;
114  }
115  else
116  _new_weights[qp] = 0.0;
117  }
118  if (!nonzero)
119  {
120  // Set the weights to a small value (1e-3) to avoid having DOFs
121  // with zeros on the diagonal, which occurs for nodes that are
122  // connected to no physical material.
123  for (unsigned qp = 0; qp < qrule->n_points(); ++qp)
124  _new_weights[qp] = 1e-3;
125  }
126  break;
127  }
128  default:
129  mooseError("Undefined option for XFEM_QRULE");
130  }
131  _have_weights = true;
132 }
133 
134 bool
135 XFEMCutElem::isPointPhysical(const Point & p) const
136 {
137  // determine whether the point is inside the physical domain of a partial element
138  bool physical_flag = true;
139  unsigned int n_cut_planes = numCutPlanes();
140  for (unsigned int plane_id = 0; plane_id < n_cut_planes; ++plane_id)
141  {
142  Point origin = getCutPlaneOrigin(plane_id);
143  Point normal = getCutPlaneNormal(plane_id);
144  Point origin2qp = p - origin;
145  if (origin2qp * normal > 0.0)
146  {
147  physical_flag = false; // Point outside pysical domain
148  break;
149  }
150  }
151  return physical_flag;
152 }
XFEM_QRULE
Definition: XFEM.h:36
unsigned int _n_qpoints
Definition: XFEMCutElem.h:36
virtual ~XFEMCutElem()
Definition: XFEMCutElem.C:28
std::vector< Point > _qp_points
Definition: XFEMCutElem.h:38
Real _physical_volfrac
Definition: XFEMCutElem.h:41
virtual Point getCutPlaneOrigin(unsigned int plane_id, MeshBase *displaced_mesh=NULL) const =0
Real getPhysicalVolumeFraction() const
Definition: XFEMCutElem.C:31
virtual Point getCutPlaneNormal(unsigned int plane_id, MeshBase *displaced_mesh=NULL) const =0
bool isPointPhysical(const Point &p) const
Definition: XFEMCutElem.C:135
std::vector< Real > _qp_weights
Definition: XFEMCutElem.h:39
std::vector< Real > _new_weights
Definition: XFEMCutElem.h:43
bool _have_weights
Definition: XFEMCutElem.h:42
virtual unsigned int numCutPlanes() const =0
virtual void computePhysicalVolumeFraction()=0
void computeXFEMWeights(QBase *qrule, Xfem::XFEM_QRULE xfem_qrule, const MooseArray< Point > &q_points)
Definition: XFEMCutElem.C:51
unsigned int _n_nodes
Definition: XFEMCutElem.h:35
std::vector< Node * > _nodes
Definition: XFEMCutElem.h:37
Real _elem_volume
Definition: XFEMCutElem.h:40
void getWeightMultipliers(MooseArray< Real > &weights, QBase *qrule, Xfem::XFEM_QRULE xfem_qrule, const MooseArray< Point > &q_points)
Definition: XFEMCutElem.C:37
XFEMCutElem(Elem *elem, unsigned int n_qpoints)
Definition: XFEMCutElem.C:20
virtual void computeMomentFittingWeights()=0