Initial Boundary Penetration Detection

r2099

framework/include/auxkernels/PenetrationAux.h

@@ -0,0 +1,35 @@
+#ifndef PENETRATIONAUX_H
+#define PENETRATIONAUX_H
+
+#include "AuxKernel.h"
+#include "PenetrationLocator.h"
+
+
+//Forward Declarations
+class PenetrationAux;
+
+template<>
+InputParameters validParams<PenetrationAux>();
+
+/** 
+ * Constant auxiliary value
+ */
+class PenetrationAux : public AuxKernel
+{
+public:
+
+  /**
+   * Factory constructor, takes parameters so that all derived classes can be built using the same
+   * constructor.
+   */
+  PenetrationAux(std::string name, MooseSystem & moose_system, InputParameters parameters);
+
+  virtual ~PenetrationAux() {}
+
+protected:
+  virtual Real computeValue();
+
+  PenetrationLocator _penetration_locator;
+};
+
+#endif //PENETRATIONAUX_H

framework/include/utils/PenetrationLocator.h

@@ -0,0 +1,32 @@
+#ifndef PENETRATIONLOCATOR_H
+#define PENETRATIONLOCATOR_H
+
+#include "libmesh_common.h"
+
+#include <vector>
+#include <map>
+
+#include "mesh.h"
+
+class PenetrationLocator
+{
+public:
+
+  PenetrationLocator(Mesh & mesh, short int master, short int slave);
+  void detectPenetration();
+
+  Real penetrationDistance(unsigned int node_id) const;
+
+private:
+
+  Real normDistance(const Elem & side, const Point & p0);
+
+  Mesh & _mesh;
+  short int _master_boundary;
+  short int _slave_boundary;
+
+  std::map<unsigned int, std::pair<unsigned int, Real> > _penetrated_elems;
+};
+
+
+#endif //PENETRATIONLOCATOR_H

framework/src/auxkernels/ConstantAux.C

@@ -4,7 +4,7 @@ template<>
 InputParameters validParams<ConstantAux>()
 {
   InputParameters params = validParams<AuxKernel>();
-  params.set<Real>("value")=0.0;
+  params.addParam<Real>("value", 0.0, "Some constant value that can be read from the input file");
   return params;
 }
 

framework/src/auxkernels/PenetrationAux.C

@@ -0,0 +1,27 @@
+#include "PenetrationAux.h"
+#include "MooseSystem.h"
+
+#include "mesh.h"
+
+template<>
+InputParameters validParams<PenetrationAux>()
+{
+  InputParameters params = validParams<AuxKernel>();
+  params.addRequiredParam<unsigned int>("master", "The boundary which will penetrate");
+  params.addRequiredParam<unsigned int>("slave", "The boundary to be penetrated");
+  return params;
+}
+
+PenetrationAux::PenetrationAux(std::string name, MooseSystem & moose_system, InputParameters parameters)
+  :AuxKernel(name, moose_system, parameters),
+   _penetration_locator(*moose_system.getMesh(), parameters.get<unsigned int>("master"), parameters.get<unsigned int>("slave"))
+{
+  // For now we are working with meshes that do not move which means we can detect penetration once!
+  _penetration_locator.detectPenetration();
+}
+
+Real
+PenetrationAux::computeValue()
+{
+  return _penetration_locator.penetrationDistance(_current_node->id());
+}

framework/src/materials/MaterialWarehouse.C

@@ -27,20 +27,28 @@ MaterialWarehouse::~MaterialWarehouse()
 std::vector<Material *>
 MaterialWarehouse::getMaterials(THREAD_ID tid, unsigned int block_id)
 {
+  std::stringstream oss;
+
   MaterialIterator mat_iter = _active_materials[tid].find(block_id);
   if (mat_iter == _active_materials[tid].end())
-    mooseError("Active Material Missing\n");
-
+  {
+    oss << "Active Material Missing for block: " << block_id << "\n";
+    mooseError(oss.str());
+  }
   return mat_iter->second;
 }
 
 std::vector<Material *>
 MaterialWarehouse::getBoundaryMaterials(THREAD_ID tid, unsigned int boundary_id)
 {
+  std::stringstream oss;
+
   MaterialIterator mat_iter = _active_boundary_materials[tid].find(boundary_id);
   if (mat_iter == _active_boundary_materials[tid].end())
-    mooseError("Active Boundary Material Missing\n");
-
+  {
+    oss << "Active Boundary Material Missing for boundary: " << boundary_id << "\n";
+    mooseError(oss.str());
+  }
   return mat_iter->second;
 }
 

framework/src/utils/PenetrationLocator.C

@@ -0,0 +1,134 @@
+#include "PenetrationLocator.h"
+
+#include "boundary_info.h"
+#include "elem.h"
+#include "plane.h"
+
+PenetrationLocator::PenetrationLocator(Mesh & mesh, short int master, short int slave)
+  : _mesh(mesh),
+    _master_boundary(master),
+    _slave_boundary(slave)
+{}
+
+
+void
+PenetrationLocator::detectPenetration()
+{
+  // Data structures to hold the element boundary information
+  std::vector< unsigned int > elem_list;
+  std::vector< unsigned short int > side_list;
+  std::vector< short int > id_list;
+
+  // Retrieve the Element Boundary data structures from the mesh
+  _mesh.boundary_info->build_side_list(elem_list, side_list, id_list);
+
+  // Data strcutres to hold the Nodal Boundary conditions
+  std::vector< unsigned int > node_list;
+  std::vector< short int > node_boundary_list;
+  _mesh.boundary_info->build_node_list_from_side_list();
+  _mesh.boundary_info->build_node_list(node_list, node_boundary_list);
+
+  // Iterator pairs for use in range functions
+  std::pair<std::vector<short int>::iterator, std::vector<short int>::iterator> bounds_node_boundary;
+  std::pair<std::vector<unsigned int>::iterator, std::vector<unsigned int>::iterator> bounds_nodes, bounds_elems;
+  std::pair<std::vector<unsigned short int>::iterator, std::vector<unsigned short int>::iterator> bounds_sides;
+
+  // Find the iterators for the master sideset nodes
+  bounds_node_boundary = std::equal_range(node_boundary_list.begin(), node_boundary_list.end(), _master_boundary);
+
+  // Bound the vector returned from the boundary_info class where the corresponding id vector matches the
+  // master boundary.  We are going to see if any of the nodes on this boundary have penetrated through
+  // elements on the slave boundary
+  bounds_nodes = std::make_pair(
+    node_list.begin() + int(bounds_node_boundary.first - node_boundary_list.begin()),
+    node_list.begin() + int(bounds_node_boundary.second - node_boundary_list.begin()));
+
+  MeshBase::const_node_iterator node = _mesh.local_nodes_begin();
+  MeshBase::const_node_iterator end_node = _mesh.local_nodes_end();
+  for ( ; node != end_node ; ++node)
+  {
+    // This is a node on the master boundary
+    if (std::binary_search(bounds_nodes.first, bounds_nodes.second, (*node)->id()))
+    {
+      // Now we need to loop over the active elements to see if we have penetrated any
+      // TODO: This might need some cleanup in parallel
+      MeshBase::const_element_iterator el = _mesh.active_local_elements_begin();
+      MeshBase::const_element_iterator end_el = _mesh.active_local_elements_end();
+      for ( ; el != end_el ; ++el)
+      {
+        const Elem* elem = *el;
+
+        bounds_node_boundary = std::equal_range(id_list.begin(), id_list.end(), _slave_boundary);
+
+        bounds_elems = std::make_pair(
+          elem_list.begin() + int(bounds_node_boundary.first - id_list.begin()),
+          elem_list.begin() + int(bounds_node_boundary.second - id_list.begin()));
+
+        std::vector<unsigned int>::iterator pos = std::lower_bound (bounds_elems.first,
+                                                                    bounds_elems.second,
+                                                                    elem->id());
+
+        // Found the local element on the slave boundary
+        if (pos != bounds_elems.second && *pos == elem->id())
+        {
+          if (elem->contains_point(**node))
+          {
+            unsigned int side_num = *(side_list.begin() + int(pos - elem_list.begin()));
+
+#ifdef DEBUG            
+            std::cout << "Node " << (*node)->id() << " contained in " << elem->id()
+                      << " through side " << side_num
+                      << ". Distance: " << normDistance(*(elem->build_side(side_num)), **node) << "\n";
+#endif            
+
+            _penetrated_elems[(*node)->id()] = std::make_pair(elem->id(), normDistance(*(elem->build_side(side_num)), **node));
+
+          }
+        }
+      }
+    }
+  }
+}
+
+Real
+PenetrationLocator::penetrationDistance(unsigned int node_id) const
+{
+  std::map<unsigned int, std::pair<unsigned int, Real> >::const_iterator found_it;
+
+  if ((found_it = _penetrated_elems.find(node_id)) != _penetrated_elems.end())
+    return found_it->second.second;
+  else
+    return 0;
+}
+
+Real
+PenetrationLocator::normDistance(const Elem & side, const Point & p0)
+{
+  Real d;
+  unsigned int dim = _mesh.mesh_dimension();
+
+  if (dim == 2)
+  {
+//    libmesh_assert(side->n_points() == 2);
+
+    Point p1 = side.point(0);
+    Point p2 = side.point(1);
+
+//    std::cerr << "\nLine Segment: (" << p1(0) << "," << p1(1) << ") - (" << p2(0) << "," << p2(1) << ") "
+//ç              << "Point: (" << p0(0) << "," << p0(1) << ")\n";
+
+
+    d = std::sqrt(std::pow(((p2(1)-p1(1))*(p0(0)-p1(0))-(p2(0)-p1(0))*(p0(1)-p1(1))),2) /
+                   (std::pow(p2(0)-p1(0),2) + std::pow(p2(1)-p1(1),2)));
+  }
+  else if (dim == 3)
+  {
+//    libmesh_assert(side.n_points() >= 3);
+
+    Plane p = Plane(side.point(0), side.point(1), side.point(2));
+
+    d = (p0 - p.closest_point(p0)).size();
+  }
+
+  return d;
+}

test/src/base/MooseTest.C

@@ -19,6 +19,7 @@
 #include "PolyReaction.h" //including our polynomial Convection Kernel
 #include "PolyForcing.h" //including our polynomial Reaction Kernel
 #include "PolyConstantAux.h" //including our polynomial Aux Kernel
+#include "PenetrationAux.h"
 //Local Includes
 #include "DiffMKernel.h"
 
@@ -59,6 +60,9 @@ namespace MooseTest
     //Registering our Aux Kernel
     AuxFactory::instance()->registerAux<MMSConstantAux>("MMSConstantAux");
 
-    AuxFactory::instance()->registerAux<PolyConstantAux>("PolyConstantAux");}
+    AuxFactory::instance()->registerAux<PolyConstantAux>("PolyConstantAux");
+
+    AuxFactory::instance()->registerAux<PenetrationAux>("PenetrationAux");
+  }
 }
 

test/tests/penetration_locator_test/penetration_locator_test.i

@@ -0,0 +1,97 @@
+[Mesh]
+  dim = 2
+  file = 2dcontact_collide.e
+#  uniform_refine = 1
+[]
+
+[Variables]
+  active = 'u'
+
+  [./u]
+    order = FIRST
+    family = LAGRANGE
+  [../]
+[]
+
+[AuxVariables]
+  [./penetration]
+    order = FIRST
+    family = LAGRANGE
+  [../]
+[]
+
+[Kernels]
+  active = 'diff'
+
+  [./diff]
+    type = Diffusion
+    variable = u
+  [../]
+[]
+
+[AuxKernels]
+  active = 'penetrate'
+
+  [./penetrate]
+    type = PenetrationAux
+    variable = penetration
+    master = 2
+    slave = 3
+  [../]
+[]
+
+[BCs]
+  active = 'block1_left block1_right block2_left block2_right'
+
+  [./block1_left]
+    type = DirichletBC
+    variable = u
+    boundary = 1
+    value = 0
+  [../]
+
+  [./block1_right]
+    type = DirichletBC
+    variable = u
+    boundary = 2
+    value = 1
+  [../]
+
+  [./block2_left]
+    type = DirichletBC
+    variable = u
+    boundary = 3
+    value = 0
+  [../]
+
+  [./block2_right]
+    type = DirichletBC
+    variable = u
+    boundary = 4
+    value = 1
+  [../]
+[]
+
+[Materials]
+  active = empty
+
+  [./empty]
+    type = EmptyMaterial
+    block = '1 2'
+  [../]
+[]
+
+[Executioner]
+  type = Steady
+  perf_log = true
+  petsc_options = '-snes_mf_operator'
+[]
+
+[Output]
+  file_base = out
+  output_initial = true
+  interval = 1
+  exodus = true
+[]
+   
+    

test/tests/penetration_locator_test/penetration_locator_test.py

@@ -0,0 +1,4 @@
+import tools
+
+def test(dofs=0, np=0):
+  tools.executeAppAndDiff(__file__,'penetration_locator_test.i',['out.e'], dofs, np)