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 "XFEMCutElem3D.h"
      11             : 
      12             : #include "EFANode.h"
      13             : #include "EFAFace.h"
      14             : #include "EFAFragment3D.h"
      15             : #include "EFAFuncs.h"
      16             : #include "XFEMFuncs.h"
      17             : #include "MooseError.h"
      18             : 
      19             : #include "libmesh/mesh.h"
      20             : #include "libmesh/elem.h"
      21             : #include "libmesh/node.h"
      22             : 
      23        2680 : XFEMCutElem3D::XFEMCutElem3D(Elem * elem,
      24             :                              const EFAElement3D * const CEMelem,
      25             :                              unsigned int n_qpoints,
      26        2680 :                              unsigned int n_sides)
      27        2680 :   : XFEMCutElem(elem, n_qpoints, n_sides), _efa_elem3d(CEMelem, true)
      28             : {
      29        2680 :   computePhysicalVolumeFraction();
      30        2680 :   computeMomentFittingWeights();
      31        2680 : }
      32             : 
      33        5360 : XFEMCutElem3D::~XFEMCutElem3D() {}
      34             : 
      35             : Point
      36     3107434 : XFEMCutElem3D::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     3107434 :   _efa_elem3d.getMasterInfo(CEMnode, master_nodes, master_weights);
      44     9164622 :   for (unsigned int i = 0; i < master_nodes.size(); ++i)
      45             :   {
      46     6057188 :     if (master_nodes[i]->category() == EFANode::N_CATEGORY_LOCAL_INDEX)
      47             :     {
      48     6057188 :       Node * node = _nodes[master_nodes[i]->id()];
      49     6057188 :       if (displaced_mesh)
      50     3595168 :         node = displaced_mesh->node_ptr(node->id());
      51     6057188 :       Point node_p((*node)(0), (*node)(1), (*node)(2));
      52     6057188 :       master_points.push_back(node_p);
      53             :     }
      54             :     else
      55           0 :       mooseError("master nodes must be local");
      56             :   }
      57     9164622 :   for (unsigned int i = 0; i < master_nodes.size(); ++i)
      58             :     node_coor += master_weights[i] * master_points[i];
      59     3107434 :   return node_coor;
      60     3107434 : }
      61             : 
      62             : void
      63       58208 : XFEMCutElem3D::computePhysicalVolumeFraction()
      64             : {
      65             :   Real frag_vol = 0.0;
      66             : 
      67             :   // collect fragment info needed by polyhedron_volume_3d()
      68             :   std::vector<std::vector<unsigned int>> frag_face_indices;
      69             :   std::vector<EFANode *> frag_nodes;
      70       58208 :   _efa_elem3d.getFragment(0)->getNodeInfo(frag_face_indices, frag_nodes);
      71       58208 :   int face_num = frag_face_indices.size();
      72       58208 :   int node_num = frag_nodes.size();
      73             : 
      74             :   int order_max = 0;
      75       58208 :   int * order = new int[face_num];
      76      355021 :   for (int i = 0; i < face_num; ++i)
      77             :   {
      78      296813 :     if (frag_face_indices[i].size() > (unsigned int)order_max)
      79       63760 :       order_max = frag_face_indices[i].size();
      80      296813 :     order[i] = frag_face_indices[i].size();
      81             :   }
      82             : 
      83       58208 :   double * coord = new double[3 * node_num];
      84      418994 :   for (unsigned int i = 0; i < frag_nodes.size(); ++i)
      85             :   {
      86      360786 :     Point p = getNodeCoordinates(frag_nodes[i]);
      87      360786 :     coord[3 * i + 0] = p(0);
      88      360786 :     coord[3 * i + 1] = p(1);
      89      360786 :     coord[3 * i + 2] = p(2);
      90             :   }
      91             : 
      92       58208 :   int * node = new int[face_num * order_max];
      93       58208 :   Xfem::i4vec_zero(face_num * order_max, node);
      94      355021 :   for (int i = 0; i < face_num; ++i)
      95     1379179 :     for (unsigned int j = 0; j < frag_face_indices[i].size(); ++j)
      96     1082366 :       node[order_max * i + j] = frag_face_indices[i][j];
      97             : 
      98             :   // compute fragment volume and volume fraction
      99       58208 :   frag_vol = Xfem::polyhedron_volume_3d(coord, order_max, face_num, node, node_num, order);
     100       58208 :   _physical_volfrac = frag_vol / _elem_volume;
     101             : 
     102       58208 :   delete[] order;
     103       58208 :   delete[] coord;
     104       58208 :   delete[] node;
     105       58208 : }
     106             : 
     107             : void
     108          60 : XFEMCutElem3D::computePhysicalFaceAreaFraction(unsigned int side)
     109             : {
     110             :   Real frag_surf = 0.0;
     111             : 
     112             :   std::vector<std::vector<unsigned int>> frag_face_indices;
     113             :   std::vector<EFANode *> frag_nodes;
     114          60 :   _efa_elem3d.getFragment(0)->getNodeInfo(frag_face_indices, frag_nodes);
     115          60 :   int face_num = frag_face_indices.size();
     116             : 
     117          60 :   EFAFace * efa_face = _efa_elem3d.getFace(side);
     118             :   bool contains_all = true;
     119             : 
     120             :   /// find a fragment surface which is covered by element side
     121         180 :   for (int i = 0; i < face_num; ++i)
     122             :   {
     123             :     contains_all = true;
     124         420 :     for (unsigned int j = 0; j < frag_face_indices[i].size(); ++j)
     125             :     {
     126         360 :       EFANode * efa_node = frag_nodes[frag_face_indices[i][j]];
     127         360 :       if (!efa_face->containsNode(efa_node))
     128             :       {
     129             :         contains_all = false;
     130             :         break;
     131             :       }
     132             :     }
     133             : 
     134         180 :     if (contains_all)
     135             :     {
     136         300 :       for (unsigned int j = 0; j < frag_face_indices[i].size(); ++j)
     137             :       {
     138         240 :         unsigned int m = ((j + 1) == frag_face_indices[i].size()) ? 0 : j + 1;
     139         240 :         Point edge_p1 = getNodeCoordinates(frag_nodes[frag_face_indices[i][j]]);
     140         240 :         Point edge_p2 = getNodeCoordinates(frag_nodes[frag_face_indices[i][m]]);
     141             : 
     142         240 :         frag_surf += 0.5 * (edge_p1(0) - edge_p2(0)) * (edge_p1(1) + edge_p2(1));
     143             :       }
     144          60 :       _physical_areafrac[side] = std::abs(frag_surf) / _elem_side_area[side];
     145             :       return;
     146             :     }
     147             :   }
     148           0 :   _physical_areafrac[side] = 1.0;
     149          60 : }
     150             : 
     151             : void
     152        2680 : XFEMCutElem3D::computeMomentFittingWeights()
     153             : {
     154             :   // TODO: 3D moment fitting method nod coded yet - use volume fraction for now
     155        2680 :   _new_weights.resize(_n_qpoints, _physical_volfrac);
     156        2680 : }
     157             : 
     158             : Point
     159      315750 : XFEMCutElem3D::getCutPlaneOrigin(unsigned int plane_id, MeshBase * displaced_mesh) const
     160             : {
     161             :   Point orig(0.0, 0.0, 0.0);
     162             :   std::vector<std::vector<EFANode *>> cut_plane_nodes;
     163     1912832 :   for (unsigned int i = 0; i < _efa_elem3d.getFragment(0)->numFaces(); ++i)
     164             :   {
     165     1597082 :     if (_efa_elem3d.getFragment(0)->isFaceInterior(i))
     166             :     {
     167      315750 :       EFAFace * face = _efa_elem3d.getFragment(0)->getFace(i);
     168             :       std::vector<EFANode *> node_line;
     169     1414206 :       for (unsigned int j = 0; j < face->numNodes(); ++j)
     170     1098456 :         node_line.push_back(face->getNode(j));
     171      315750 :       cut_plane_nodes.push_back(node_line);
     172      315750 :     }
     173             :   }
     174      315750 :   if (cut_plane_nodes.size() == 0)
     175           0 :     mooseError("no cut plane found in this element");
     176      315750 :   if (plane_id < cut_plane_nodes.size()) // valid plane_id
     177             :   {
     178             :     std::vector<Point> cut_plane_points;
     179      720078 :     for (unsigned int i = 0; i < cut_plane_nodes[plane_id].size(); ++i)
     180     1119216 :       cut_plane_points.push_back(getNodeCoordinates(cut_plane_nodes[plane_id][i], displaced_mesh));
     181             : 
     182      720078 :     for (unsigned int i = 0; i < cut_plane_points.size(); ++i)
     183             :       orig += cut_plane_points[i];
     184      160470 :     orig /= cut_plane_points.size();
     185      160470 :   }
     186      315750 :   return orig;
     187      315750 : }
     188             : 
     189             : Point
     190      519468 : XFEMCutElem3D::getCutPlaneNormal(unsigned int plane_id, MeshBase * displaced_mesh) const
     191             : {
     192             :   Point normal(0.0, 0.0, 0.0);
     193             :   std::vector<std::vector<EFANode *>> cut_plane_nodes;
     194     3336602 :   for (unsigned int i = 0; i < _efa_elem3d.getFragment(0)->numFaces(); ++i)
     195             :   {
     196     2817134 :     if (_efa_elem3d.getFragment(0)->isFaceInterior(i))
     197             :     {
     198      519468 :       EFAFace * face = _efa_elem3d.getFragment(0)->getFace(i);
     199             :       std::vector<EFANode *> node_line;
     200     2431400 :       for (unsigned int j = 0; j < face->numNodes(); ++j)
     201     1911932 :         node_line.push_back(face->getNode(j));
     202      519468 :       cut_plane_nodes.push_back(node_line);
     203      519468 :     }
     204             :   }
     205      519468 :   if (cut_plane_nodes.size() == 0)
     206           0 :     mooseError("no cut plane found in this element");
     207      519468 :   if (plane_id < cut_plane_nodes.size()) // valid plane_id
     208             :   {
     209             :     std::vector<Point> cut_plane_points;
     210     1737272 :     for (unsigned int i = 0; i < cut_plane_nodes[plane_id].size(); ++i)
     211     2746168 :       cut_plane_points.push_back(getNodeCoordinates(cut_plane_nodes[plane_id][i], displaced_mesh));
     212             : 
     213             :     Point center(0.0, 0.0, 0.0);
     214     1737272 :     for (unsigned int i = 0; i < cut_plane_points.size(); ++i)
     215             :       center += cut_plane_points[i];
     216      364188 :     center /= cut_plane_points.size();
     217             : 
     218     1737272 :     for (unsigned int i = 0; i < cut_plane_points.size(); ++i)
     219             :     {
     220     1373084 :       unsigned int iplus1 = i < cut_plane_points.size() - 1 ? i + 1 : 0;
     221             :       Point ray1 = cut_plane_points[i] - center;
     222     1373084 :       Point ray2 = cut_plane_points[iplus1] - center;
     223     1373084 :       normal += ray1.cross(ray2);
     224             :     }
     225             :     normal /= cut_plane_points.size();
     226      364188 :   }
     227      519468 :   Xfem::normalizePoint(normal);
     228      519468 :   return normal;
     229      519468 : }
     230             : 
     231             : void
     232           0 : XFEMCutElem3D::getCrackTipOriginAndDirection(unsigned /*tip_id*/,
     233             :                                              Point & /*origin*/,
     234             :                                              Point & /*direction*/) const
     235             : {
     236             :   // TODO: not implemented for 3D
     237           0 :   mooseError("getCrackTipOriginAndDirection not yet implemented for XFEMCutElem3D");
     238             : }
     239             : 
     240             : void
     241           0 : XFEMCutElem3D::getFragmentFaces(std::vector<std::vector<Point>> & /*frag_faces*/,
     242             :                                 MeshBase * /*displaced_mesh*/) const
     243             : {
     244             :   // TODO: not implemented for 3D
     245           0 :   mooseError("getFragmentFaces not yet implemented for XFEMCutElem3D");
     246             : }
     247             : 
     248             : const EFAElement *
     249      930434 : XFEMCutElem3D::getEFAElement() const
     250             : {
     251      930434 :   return &_efa_elem3d;
     252             : }
     253             : 
     254             : unsigned int
     255        5190 : XFEMCutElem3D::numCutPlanes() const
     256             : {
     257             :   unsigned int counter = 0;
     258       36416 :   for (unsigned int i = 0; i < _efa_elem3d.getFragment(0)->numFaces(); ++i)
     259       31226 :     if (_efa_elem3d.getFragment(0)->isFaceInterior(i))
     260        5190 :       counter += 1;
     261        5190 :   return counter;
     262             : }
     263             : 
     264             : void
     265      203718 : XFEMCutElem3D::getIntersectionInfo(unsigned int plane_id,
     266             :                                    Point & normal,
     267             :                                    std::vector<Point> & intersectionPoints,
     268             :                                    MeshBase * displaced_mesh) const
     269             : {
     270      203718 :   intersectionPoints.clear();
     271             :   std::vector<std::vector<EFANode *>> cut_plane_nodes;
     272     1423770 :   for (unsigned int i = 0; i < _efa_elem3d.getFragment(0)->numFaces(); ++i)
     273             :   {
     274     1220052 :     if (_efa_elem3d.getFragment(0)->isFaceInterior(i))
     275             :     {
     276      203718 :       EFAFace * face = _efa_elem3d.getFragment(0)->getFace(i);
     277             :       std::vector<EFANode *> node_line;
     278     1017194 :       for (unsigned int j = 0; j < face->numNodes(); ++j)
     279      813476 :         node_line.push_back(face->getNode(j));
     280      203718 :       cut_plane_nodes.push_back(node_line);
     281      203718 :     }
     282             :   }
     283      203718 :   if (cut_plane_nodes.size() == 0)
     284           0 :     mooseError("No cut plane found in this element");
     285             : 
     286      203718 :   if (plane_id < cut_plane_nodes.size()) // valid plane_id
     287             :   {
     288      203718 :     intersectionPoints.resize(cut_plane_nodes[plane_id].size());
     289     1017194 :     for (unsigned int i = 0; i < cut_plane_nodes[plane_id].size(); ++i)
     290      813476 :       intersectionPoints[i] = getNodeCoordinates(cut_plane_nodes[plane_id][i], displaced_mesh);
     291             :   }
     292             : 
     293      203718 :   normal = getCutPlaneNormal(plane_id, displaced_mesh);
     294      203718 : }