mortaring_helpers.h

Go to the documentation of this file.

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
#ifndef MORTARING_HELPERS_HEADER
#define MORTARING_HELPERS_HEADER
 
#include "mesh.h"
#include "Oomph/Constraints/constraint_elements.h"
 
using namespace oomph;
 
 
// From a list of nodes which are to be mortared create a list of spatially unique nodes so that we don't create singular constraints
// use point elements to perform setup multi domain interactions 
// Do this twice, grabbing the element and local coordiante for each point element for each mesh
// Build constraint elements and return them
 
// This doesn't use a mortaring status/dof - although this could be used for other implementations I suppose but it can be separate to this process
 
// setup_mortaring_constraints_at_nodes attempt to match at nodes first before doing locate zetas.
// for any nodes which haven't been able to be identified as having a partner node try to locate zetas
 
namespace MortaringHelpers
{
 
// A helper class which merges two meshes and computes boundary elements on the new boundary numbering scheme
// Not used anywhere.
class MergedSolidMesh : public virtual SolidMesh
{
public:
 MergedSolidMesh(Vector<SolidMesh*> solid_mesh_pt) : SolidMesh(solid_mesh_pt)
 {
  // Solid mesh constructor has merged meshes
 
  // We now need to update our representation of the boundary numbers
 
  // Loop over nodes adding them to the appropriate boundary lookup schemes
  // in the mesh. Shift by the required number of boundaries
  const unsigned n_mesh = solid_mesh_pt.size();
  Vector<unsigned> Boundary_shift(n_mesh);
 
  unsigned n_boundary = 0;
  for(unsigned m=0; m<n_mesh; m++)
  {
   n_boundary += solid_mesh_pt[m]->nboundary();
   if(m==n_mesh-1) continue;
   Boundary_shift[m+1] = solid_mesh_pt[m]->nboundary();
  }
  this->set_nboundary(n_boundary);
  Boundary_element_pt.resize(n_boundary);
  Face_index_at_boundary.resize(n_boundary);
 
  for(unsigned m=0; m<n_mesh; m++)
  {
   const unsigned n_node = solid_mesh_pt[m]->nnode();
   for (unsigned nd = 0; nd < n_node; nd++)
   {
    Node* nd_pt = solid_mesh_pt[m]->node_pt(nd);
    if (nd_pt->is_on_boundary())
    {
     // Get set of boundaries that the node is on
     std::set<unsigned>* boundaries_pt;
     nd_pt->get_boundaries_pt(boundaries_pt);
 
     // Add node pointer to the appropriate Boundary_node_pt vectors
     std::set<unsigned>::const_iterator it;
     for (it = boundaries_pt->begin(); it != boundaries_pt->end(); it++)
     {
      Boundary_node_pt[*it + Boundary_shift[m]].push_back(nd_pt);
     }
    }
   }
 
   // Get the list of elements on the boundary of the constituent mesh and add them to the new shifted boundary
   for(unsigned b=0; b<solid_mesh_pt[m]->nboundary(); b++)
   {
    Vector<FiniteElement*> boundary_elements(solid_mesh_pt[m]->nboundary_element(b));
    Vector<int> boundary_faces(solid_mesh_pt[m]->nboundary_element(b));
    for(unsigned e=0;e<solid_mesh_pt[m]->nboundary_element(b); e++)
    {
     boundary_elements[e] = solid_mesh_pt[m]->boundary_element_pt(b,e);
     boundary_faces[e] = solid_mesh_pt[m]->face_index_at_boundary(b,e);
    }
    Boundary_element_pt[Boundary_shift[m]+b] = boundary_elements;
    Face_index_at_boundary[Boundary_shift[m]+b] = boundary_faces;
   }
  }
 }
};
 
class PointIntegralScheme : public virtual Integral
{
public:
 PointIntegralScheme(){};
 
 PointIntegralScheme(const PointIntegralScheme& dummy) = delete;
 
 void operator=(const PointIntegralScheme&) = delete;
 
 unsigned nweight() const
 {
  return 1;
 }
 
 double knot(const unsigned& i, const unsigned& j) const
 {
  return 0.0;
 }
 
 double weight(const unsigned& i) const
 {
  return 1.0;
 }
};
 
 
// A helper element for locating the point in a mesh at which a node in another mesh sits
class PointElementWithExternalElement : public virtual PointElement, public virtual ElementWithExternalElement
{
private:
 // Default integral scheme, defines the locations of the integral points
 static PointIntegralScheme Interface_default_integration_scheme;
 
public:
 PointElementWithExternalElement(const Vector<double>& coord, const unsigned& n_interaction=1)
 {
  // Set the dimensions of the element and the nodes
  this->set_dimension(coord.size());
  // By default there is 1 interaction
  this->set_ninteraction(n_interaction);
  // Set the integration scheme
  this->set_integration_scheme(&Interface_default_integration_scheme);
 
  // Position the node at the particle
  Node* node_pt = this->construct_node(0, dynamic_cast<TimeStepper*>(&Mesh::Default_TimeStepper));
  for(unsigned i=0; i<coord.size(); i++)
  {
   node_pt->x(i) = coord[i];
  }
 }
 
 void fill_in_contribution_to_jacobian(Vector<double>&  residuals, DenseMatrix<double>& jacobian){};
};
PointIntegralScheme PointElementWithExternalElement::Interface_default_integration_scheme;
 
 
 
 
 
 
 
 
 
 
static double Mortaring_Constraints_Distance_Tolerance = 1e-10;
// Setup mortaring constraint interactions between two meshes at the locations of the nodes
// initially looks for co-located node pairs then with any remaining nodes uses locate zeta
// to find a suitable co-locating point in another element
// Provide two vectors, one for the nodes in each element
// For efficiency assumes that only 2 nodes will be co-located and will not check nodes against nodes of the same mesh
template<class NodeNodeConstraintClass, class NodeElementConstraintClass, class ELEMENTTYPE>
Vector<ConstraintElement*> setup_constraint_elements_at_nodes(Problem* problem_pt,
                                                              Vector<Vector<Node*>> node_pt,
                                                              Vector<SolidMesh*> mesh_pt,
                                                              const bool& accept_failure = false)
{
 // Vector of constraint elements built
 Vector<ConstraintElement*> constraint_element_pt;
 
 // Assume all the nodes/elements are the same dimension
 const unsigned dim = node_pt[0][0]->ndim();
 
 // The number of nodes from each mesh
 const Vector<long unsigned> n_node = {node_pt[0].size(), node_pt[1].size()};
 
 // Allocate for nodes which were not matched to other nodes
 Vector<Vector<unsigned>> node_was_matched(2);
 node_was_matched[0].resize(n_node[0], 0);
 node_was_matched[1].resize(n_node[1], 0);
 
 // First attempt to match each node to another in the other mesh
 unsigned n_matched_nodes = 0;
 for(unsigned i=0; i<n_node[0]; i++)
 {
  SolidNode* node_0_pt = dynamic_cast<SolidNode*>(node_pt[0][i]);
  for(unsigned j=0; j<n_node[1]; j++)
  {
   // If the node has already been matched to another node
   if(node_was_matched[1][j]) continue;
 
   SolidNode* node_1_pt = dynamic_cast<SolidNode*>(node_pt[1][j]);
 
   // get the distance between the nodes TODO is it faster to just check absolute distance in each spatial direction?
   double distance = 0.0;
   for(unsigned d=0; d<dim; d++)
   {
    distance += (node_0_pt->xi(d) - node_1_pt->xi(d))*
                (node_0_pt->xi(d) - node_1_pt->xi(d));
   }
   // if co-located build constraint element
   if(sqrt(distance) < Mortaring_Constraints_Distance_Tolerance)
   {
    // If the nodes are the same then skip this, we can't mortar a node to itself
    // This check is useful because it will allow a mesh to be mortared to itself
    if(node_0_pt == node_1_pt) continue;
    // Record that these nodes were matched
    node_was_matched[0][i] = 1;
    node_was_matched[1][j] = 1;
    n_matched_nodes++;
    // Build the node-node constraint element
    constraint_element_pt.push_back(new NodeNodeConstraintClass(node_0_pt, node_1_pt));
 
    // stop checking for this node in the first mesh
    break;
   }
  }
 }
 
 // If we have matched all nodes then return
 if((n_node[0] - n_matched_nodes) == 0 && (n_node[1] - n_matched_nodes) == 0 ) return constraint_element_pt;
 
 // If we get to here then we know that some nodes have not been mortared.
 // If the mesh pointers are the same then we know we have failed. We can't perform the following procedure
 // if the meshes are the same.
 if(mesh_pt[0] == mesh_pt[1])
 {
  throw OomphLibError("Mortaring failed - meshes are the same and node-node mortaring has not captured all nodes",
                   OOMPH_CURRENT_FUNCTION,
                   OOMPH_EXCEPTION_LOCATION);
 }
 
 const bool Previous_Accept_failed_locate_zeta_in_setup_multi_domain_interaction
  = Multi_domain_functions::Accept_failed_locate_zeta_in_setup_multi_domain_interaction;
 Multi_domain_functions::Accept_failed_locate_zeta_in_setup_multi_domain_interaction = accept_failure;
 
                                                                       
 // Figure out how many nodes are left to be matched from each mesh
 Vector<long unsigned> n_unmatched_nodes{n_node[0] - n_matched_nodes,
                                         n_node[1] - n_matched_nodes};
 
 for(unsigned m : {0,1})
 {
  if(n_unmatched_nodes[m] == 0) continue;
  
  // Get the nodes which were not matched
  Vector<Node*> unmatched_node_pt(n_unmatched_nodes[m]);
  unsigned count = 0;
  for(unsigned i=0; i<n_node[m]; i++)
  {
   if(!node_was_matched[m][i])
   {
    unmatched_node_pt[count++] = node_pt[m][i];
   }
  }
 
  // Create a mesh for locating the zetas in the other mesh
  Mesh* Locating_mesh_pt = new Mesh;
 
  Vector<PointElementWithExternalElement*> locating_elements_pt(n_unmatched_nodes[m]);
  for(unsigned i=0; i<n_unmatched_nodes[m]; i++)
  {
   Vector<double> x(dim);
   for(unsigned d=0; d<dim; d++)
   {
    x[d] = dynamic_cast<SolidNode*>(unmatched_node_pt[i])->xi(d);
   }
 
   locating_elements_pt[i] = new PointElementWithExternalElement(x);
   Locating_mesh_pt->add_element_pt(locating_elements_pt[i]);
  }
 
  // Setup multi domain interactions with the second mesh
  Multi_domain_functions::setup_multi_domain_interaction<ELEMENTTYPE>(problem_pt,
                                                                      Locating_mesh_pt,
                                                                      mesh_pt[unsigned(1-m)],
                                                                      0);
 
  // Loop over the locating elements
  for(unsigned i=0; i<n_unmatched_nodes[m]; i++)
  {
   // If the external element was not found skip this element - no need to crash out, that would have happened
   // earlier if this was an issue (if accept_failure == false)
   if(locating_elements_pt[i]->external_element_pt(0,0) == nullptr) continue;
   constraint_element_pt.push_back(new NodeElementConstraintClass(dynamic_cast<SolidNode*>(unmatched_node_pt[i]),
                                                                  locating_elements_pt[i]->external_element_pt(0,0),
                                                                  locating_elements_pt[i]->external_element_local_coord(0,0)));
  }
 }
 Multi_domain_functions::Accept_failed_locate_zeta_in_setup_multi_domain_interaction = Previous_Accept_failed_locate_zeta_in_setup_multi_domain_interaction;
 
 return constraint_element_pt;
}
 
 
 
 
 
// ////////////////////////////////////////////////////////////////////
// Required to find and remove redundant mortaring elements from a problem
// Currently only works for nodes.
// To perform the same procedure on general Node-Element mortaring elements
//  I think the Jacobian will have to be interrogated to remove linearly
//  dependent rows.
// ////////////////////////////////////////////////////////////////////
 
// Uses the depth first search (DFS)
 
// Enum to represent the colours of the nodes in the DFS (WHITE=unvisited, GREY=being visited, BLACK=fully processed)
enum Colour { WHITE, GREY, BLACK };
 
// Depth First Search (DFS) to detect cycles and identify redundant branches
template<class NodeNodeConstraintElementType>
void dfs(Node* node_pt,
   Node* parent_node_pt,
   std::map<Node*, Vector<std::pair<NodeNodeConstraintElementType*, Node*>>>& tree,
   std::map<Node*, Colour>& colours,
   Vector<NodeNodeConstraintElementType*>& pruned_branches)
{
 // Mark the current node as being visited (GREY)
 colours[node_pt] = GREY;
 
 // Get all branches (edges) from the current node
 Vector<std::pair<NodeNodeConstraintElementType*, Node*>> branches = tree[node_pt];
 for (auto& branch : branches)
 {
  Node* next_node = branch.second;
  NodeNodeConstraintElementType* element = branch.first;
 
  // If the branch has already been pruned, skip it
  if (std::find(pruned_branches.begin(), pruned_branches.end(), element) != pruned_branches.end())
  {
   continue;
  }
 
  // If we encounter a GREY node that is not the parent, a cycle is detected
  if (colours[next_node] == GREY && next_node != parent_node_pt)
  {
   pruned_branches.push_back(element);
  }
  // If the next node is unvisited (WHITE), continue DFS
  else if (colours[next_node] == WHITE)
  {
   dfs(next_node, node_pt, tree, colours, pruned_branches);
  }
 }
 
 // Mark the current node as fully processed (BLACK)
 colours[node_pt] = BLACK;
}
 
// Function to find redundant constraint elements using DFS
template<class NodeNodeConstraintElementType>
Vector<NodeNodeConstraintElementType*> find_redundant_constraint_elements(Vector<NodeNodeConstraintElementType*> element_pt)
{
 // Create a tree representation from the constraint elements
 std::map<Node*, Vector<std::pair<NodeNodeConstraintElementType*, Node*>>> tree;
 for (auto& element : element_pt)
 {
  Node* node1 = element->solid_node_pt(0);
  Node* node2 = element->solid_node_pt(1);
  tree[node1].push_back(std::make_pair(element, node2));
  tree[node2].push_back(std::make_pair(element, node1));
 }
 
 // Vector to store the pruned branches
 Vector<NodeNodeConstraintElementType*> pruned_branches;
 // Map to store the colours of the nodes
 std::map<Node*, Colour> colours;
 
 // Initialize all nodes as unvisited (WHITE)
 for (auto& node_branches : tree)
 {
  Node* root_node = node_branches.first;
  if (colours[root_node] == WHITE)
  {
   // // Perform DFS to detect cycles starting from the root node
   dfs(root_node, nullptr, tree, colours, pruned_branches);
  }
 }
 
 return pruned_branches;
}
 
 
 
}; // End namespace
 
#endif