XYCrossFieldSolver Class Reference

Solves for a cross field, the smooth field of 4-fold symmetric directions that gives XYFrontalDelaunayGenerator the local frame in which it measures the target size ahead of the front. More...

#include <XYCrossFieldSolver.h>

Public Types
using	NodalCrossField = std::map< dof_id_type, std::complex< Real > >
	Values of the complex cross field z at the nodes of the background mesh, keyed by node id.

Public Member Functions
	XYCrossFieldSolver (const libMesh::MeshBase &background_mesh, const std::map< dof_id_type, Real > &boundary_tangent_angles)
	Build a cross field solver on a background triangulation.
void	solve ()
	Assemble and solve the two Laplace problems, then normalize the nodal field and collect the singular nodes.
Real	theta (const Point &point) const
std::pair< Point, Point >	crossFrame (const Point &point) const
const NodalCrossField &	nodalCrossField () const
const std::vector< dof_id_type > &	singularNodes () const

Private Member Functions
void	assembleLaplace ()
	Assemble the Laplace operator for both variables and the nodal penalty Dirichlet conditions.
void	extractNodalCrossField ()
	Normalize the solved nodal values of z and record the nodes where z vanishes.
std::complex< Real >	interpolatedCrossField (const Point &point) const

Static Private Member Functions
static void	assembleSystem (libMesh::EquationSystems &es, const std::string &system_name)
	libMesh assembly callback; forwards to assembleLaplace() on the solver stashed in es.

Private Attributes
libMesh::Parallel::Communicator	_communicator
	Serial communicator, so that the solve is reproducible independent of the process count.
libMesh::ReplicatedMesh	_mesh
	Private copy of the background triangulation, sharing the caller's node ids.
const NodalCrossField	_boundary_cross_field
	Dirichlet value exp(4 i theta_t) of each boundary node, keyed by node id.
const libMesh::FEType	_fe_type
	Linear Lagrange, the discretization of both components of z.
std::unique_ptr< libMesh::EquationSystems >	_equation_systems
	Holds the solve; kept past solve() because shape function evaluation needs it.
libMesh::LinearImplicitSystem *	_system
	System holding both components of z; owned by _equation_systems.
unsigned int	_real_variable
	Variable number of the real part of z.
unsigned int	_imaginary_variable
	Variable number of the imaginary part of z.
std::unique_ptr< libMesh::PointLocatorBase >	_point_locator
	Locates the element containing a query point.
bool	_solved
	Whether solve() has produced the results the queries read.
NodalCrossField	_nodal_cross_field
	Unit-magnitude nodal values of z, keyed by node id.
std::vector< dof_id_type >	_singular_nodes
	Ascending ids of the nodes where z vanishes.

Detailed Description

Solves for a cross field, the smooth field of 4-fold symmetric directions that gives XYFrontalDelaunayGenerator the local frame in which it measures the target size ahead of the front.

The cross field is represented by the complex field z whose argument carries the 4-fold symmetry,

Laplace(z) = 0, z = exp(4 i theta_t) on the boundary, theta = arg(z) / 4

where theta_t is the boundary tangent angle and theta is the angle of one of the four directions of the cross at any point of the domain, which is theta_t itself on the boundary. Working with z instead of theta makes the problem linear and single valued: theta itself is only defined modulo pi/2, so it can neither be solved for nor interpolated directly. Zeros of z are the singularities of the field, the points where it has no direction; the frame is discontinuous around them, and singularNodes() reports the nodes where the solve found one.

The solve is a pair of real Laplace problems (one for the real part of z, one for the imaginary part) discretized with linear Lagrange finite elements on the background triangulation. It runs on a private serial communicator so that the field is identical on every process and for every process count.

Definition at line 53 of file XYCrossFieldSolver.h.

Member Typedef Documentation

NodalCrossField

using XYCrossFieldSolver::NodalCrossField = std::map<dof_id_type, std::complex<Real> >

Values of the complex cross field z at the nodes of the background mesh, keyed by node id.

Definition at line 57 of file XYCrossFieldSolver.h.

Constructor & Destructor Documentation

XYCrossFieldSolver()

XYCrossFieldSolver::XYCrossFieldSolver	(	const libMesh::MeshBase &	background_mesh,
		const std::map< dof_id_type, Real > &	boundary_tangent_angles
	)

Build a cross field solver on a background triangulation.

Parameters

background_mesh	Replicated background mesh of TRI3 elements covering the domain. The mesh is copied internally with its node ids preserved, so it need not outlive this object, and its node ids are the ids used by both boundary_tangent_angles and every result this class reports.
boundary_tangent_angles	Boundary tangent angle theta_t, in radians, for each boundary node of background_mesh, keyed by node id. Every node on the domain boundary must appear; nodes absent from this map are treated as interior.

Definition at line 66 of file XYCrossFieldSolver.C.

  : _communicator(MPI_COMM_SELF),
    _mesh(_communicator),
    _boundary_cross_field(boundaryCrossField(boundary_tangent_angles)),
    _fe_type(FIRST, LAGRANGE),
    _system(nullptr),
    _real_variable(0),
    _imaginary_variable(0),
    _solved(false)
{
  // Solving across a distributed mesh is not supported here: only the local portion of one would
  // be copied, which would silently solve on a fragment of the domain.
  if (!background_mesh.is_replicated())
    mooseError("XYCrossFieldSolver: The cross field background mesh must be replicated.");
 
  // Node ids key the boundary data and every reported result, so they have to survive the copy onto
  // the serial communicator.
  _mesh.allow_renumbering(false);
  _mesh.copy_nodes_and_elements(background_mesh);
  _mesh.prepare_for_use();
 
  for (const auto & elem : _mesh.element_ptr_range())
    if (elem->type() != TRI3)
      mooseError("XYCrossFieldSolver: The cross field background mesh must contain only TRI3 "
                 "elements, but it contains an element of type ",
                 Utility::enum_to_string(elem->type()),
                 ".");
 
  for (const auto & [node_id, _] : _boundary_cross_field)
    if (!_mesh.query_node_ptr(node_id))
      mooseError("XYCrossFieldSolver: A boundary tangent angle was supplied for node ",
                 node_id,
                 ", which is not a node of the cross field background mesh.");
 
  _equation_systems = std::make_unique<EquationSystems>(_mesh);
  _system = &_equation_systems->add_system<LinearImplicitSystem>("cross_field");
  _real_variable = _system->add_variable("z_real", _fe_type);
  _imaginary_variable = _system->add_variable("z_imaginary", _fe_type);
  _system->attach_assemble_function(assembleSystem);
  _equation_systems->parameters.set<XYCrossFieldSolver *>("cross_field_solver") = this;
  _equation_systems->init();
 
  _point_locator = PointLocatorBase::build(TREE_LOCAL_ELEMENTS, _mesh);
  _point_locator->enable_out_of_mesh_mode();
}

Member Function Documentation

assembleLaplace()

void XYCrossFieldSolver::assembleLaplace ( )

private

Assemble the Laplace operator for both variables and the nodal penalty Dirichlet conditions.

Definition at line 135 of file XYCrossFieldSolver.C.

{
  const DofMap & dof_map = _system->get_dof_map();
  std::unique_ptr<FEBase> fe(FEBase::build(_mesh.mesh_dimension(), _fe_type));
  QGauss qrule(_mesh.mesh_dimension(), _fe_type.default_quadrature_order());
  fe->attach_quadrature_rule(&qrule);
 
  const std::vector<Real> & JxW = fe->get_JxW();
  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
 
  DenseMatrix<Number> Ke;
  DenseVector<Number> Fe_real;
  DenseVector<Number> Fe_imaginary;
  std::vector<dof_id_type> real_dof_indices;
  std::vector<dof_id_type> imaginary_dof_indices;
 
  auto & system_matrix = _system->get_system_matrix();
 
  for (const auto & elem : _mesh.active_local_element_ptr_range())
  {
    fe->reinit(elem);
    dof_map.dof_indices(elem, real_dof_indices, _real_variable);
    dof_map.dof_indices(elem, imaginary_dof_indices, _imaginary_variable);
 
    const auto n_dofs = real_dof_indices.size();
    Ke.resize(n_dofs, n_dofs);
    Fe_real.resize(n_dofs);
    Fe_imaginary.resize(n_dofs);
 
    for (const auto qp : make_range(qrule.n_points()))
      for (const auto i : make_range(n_dofs))
        for (const auto j : make_range(n_dofs))
          Ke(i, j) += JxW[qp] * (dphi[i][qp] * dphi[j][qp]);
 
    // The Dirichlet data is given per node, so the boundary condition is a penalty on the node's
    // own equation rather than a side integral. A node shared by several elements accumulates the
    // penalty once per element on both the diagonal and the right hand side, which leaves the value
    // it pins unchanged.
    mooseAssert(n_dofs == elem->n_nodes(),
                "The nodal penalty indexes degrees of freedom by local node number, which linear "
                "Lagrange on TRI3 must supply one of per element node.");
    for (const auto i : make_range(n_dofs))
    {
      const auto boundary_value = _boundary_cross_field.find(elem->node_id(i));
      if (boundary_value == _boundary_cross_field.end())
        continue;
 
      Ke(i, i) += dirichlet_penalty;
      Fe_real(i) += dirichlet_penalty * boundary_value->second.real();
      Fe_imaginary(i) += dirichlet_penalty * boundary_value->second.imag();
    }
 
    // Both components of z obey the same Laplace operator and couple only through the boundary
    // data, so one element matrix serves both variables.
    system_matrix.add_matrix(Ke, real_dof_indices);
    system_matrix.add_matrix(Ke, imaginary_dof_indices);
    _system->rhs->add_vector(Fe_real, real_dof_indices);
    _system->rhs->add_vector(Fe_imaginary, imaginary_dof_indices);
  }
}

Referenced by assembleSystem().

assembleSystem()

void XYCrossFieldSolver::assembleSystem ( libMesh::EquationSystems &  es, const std::string &  system_name  )

staticprivate

libMesh assembly callback; forwards to assembleLaplace() on the solver stashed in es.

Definition at line 129 of file XYCrossFieldSolver.C.

{
  es.parameters.get<XYCrossFieldSolver *>("cross_field_solver")->assembleLaplace();
}

Referenced by XYCrossFieldSolver().

crossFrame()

std::pair< Point, Point > XYCrossFieldSolver::crossFrame

(

const Point &

point

)

const

Returns

The local cross frame at point as the pair (u, v), with u = (cos theta, sin theta) and v the in-plane perpendicular of u.

Definition at line 260 of file XYCrossFieldSolver.C.

{
  const Real angle = theta(point);
  const Point u(std::cos(angle), std::sin(angle), 0.0);
  return {u, Point(-u(1), u(0), 0.0)};
}

extractNodalCrossField()

void XYCrossFieldSolver::extractNodalCrossField ( )

private

Normalize the solved nodal values of z and record the nodes where z vanishes.

Definition at line 197 of file XYCrossFieldSolver.C.

{
  const auto system_number = _system->number();
  const auto & solution = *_system->solution;
 
  _nodal_cross_field.clear();
  _singular_nodes.clear();
 
  for (const auto & node : _mesh.node_ptr_range())
  {
    // The system lives on a serial communicator, so every degree of freedom is local here.
    const std::complex<Real> z(solution(node->dof_number(system_number, _real_variable, 0)),
                               solution(node->dof_number(system_number, _imaginary_variable, 0)));
    const Real magnitude = std::abs(z);
 
    if (magnitude < singular_magnitude)
      _singular_nodes.push_back(node->id());
 
    // A vanishing z carries no direction at all, so fall back on the real axis instead of dividing
    // by zero. Such a node is always singular by the test above.
    const auto normalized =
        magnitude > minimum_magnitude ? z / magnitude : std::complex<Real>(1.0, 0.0);
    _nodal_cross_field.emplace(node->id(), normalized);
  }
 
  std::sort(_singular_nodes.begin(), _singular_nodes.end());
}

Referenced by solve().

interpolatedCrossField()

std::complex< Real > XYCrossFieldSolver::interpolatedCrossField ( const Point &  point ) const

private

Returns

The cross field z interpolated over the element containing point.

Definition at line 226 of file XYCrossFieldSolver.C.

{
  mooseAssert(_solved, "solve() must be called before the cross field can be queried.");
 
  const Elem * elem = (*_point_locator)(point);
  if (!elem)
    mooseError("XYCrossFieldSolver: No element was found to contain point ", point);
 
  const Point reference_point = FEMap::inverse_map(elem->dim(), elem, point);
  FEComputeData fe_data(*_equation_systems, reference_point);
  FEInterface::compute_data(elem->dim(), _fe_type, elem, fe_data);
 
  mooseAssert(fe_data.shape.size() == elem->n_nodes(),
              "Linear Lagrange on TRI3 must supply one shape function per element node.");
 
  // Interpolating z rather than theta is what makes this well posed: theta is only defined modulo
  // pi/2, so averaging nodal angles is wrong wherever the branch cut falls between the nodes.
  std::complex<Real> interpolated(0.0, 0.0);
  for (const auto i : index_range(fe_data.shape))
  {
    const Real shape_value = fe_data.shape[i];
    interpolated += shape_value * libmesh_map_find(_nodal_cross_field, elem->node_id(i));
  }
 
  return interpolated;
}

Referenced by theta().

nodalCrossField()

const NodalCrossField & XYCrossFieldSolver::nodalCrossField ( ) const

inline

Returns

The unit-magnitude cross field value z at every node of the background mesh, keyed by node id. A node whose solved magnitude was too small to normalize is assigned z = 1; such a node is always among those reported by singularNodes().

Definition at line 97 of file XYCrossFieldSolver.h.

{ return _nodal_cross_field; }

singularNodes()

const std::vector< dof_id_type > & XYCrossFieldSolver::singularNodes ( ) const

inline

Returns

Ascending ids of the singular nodes, the nodes where the solved magnitude of z fell below 0.1 before normalization, so that the field has no direction there and the frame is discontinuous around them.

Definition at line 104 of file XYCrossFieldSolver.h.

{ return _singular_nodes; }

solve()

void XYCrossFieldSolver::solve

(

)

Assemble and solve the two Laplace problems, then normalize the nodal field and collect the singular nodes.

Must be called before any query below.

Definition at line 114 of file XYCrossFieldSolver.C.

{
  _equation_systems->parameters.set<Real>("linear solver tolerance") = linear_solver_tolerance;
  _system->solve();
 
  const auto reason = _system->get_linear_solver()->get_converged_reason();
  if (reason < 0)
    mooseError("XYCrossFieldSolver: The cross field solve failed to converge with reason: ",
               Utility::enum_to_string(reason));
 
  extractNodalCrossField();
  _solved = true;
}

theta()

Real XYCrossFieldSolver::theta

(

const Point &

point

)

const

Returns

The cross field angle theta at point, in radians, obtained by interpolating z over the containing element and taking arg(z) / 4. Since arg() returns a value in (-pi, pi], theta lands in (-pi/4, pi/4], the canonical representative of the cross direction modulo pi/2.

Definition at line 254 of file XYCrossFieldSolver.C.

{
  return std::arg(interpolatedCrossField(point)) / 4.0;
}

Referenced by crossFrame().

Member Data Documentation

_boundary_cross_field

const NodalCrossField XYCrossFieldSolver::_boundary_cross_field

private

Dirichlet value exp(4 i theta_t) of each boundary node, keyed by node id.

Definition at line 126 of file XYCrossFieldSolver.h.

Referenced by assembleLaplace(), and XYCrossFieldSolver().

_communicator

libMesh::Parallel::Communicator XYCrossFieldSolver::_communicator

private

Serial communicator, so that the solve is reproducible independent of the process count.

Definition at line 120 of file XYCrossFieldSolver.h.

_equation_systems

std::unique_ptr<libMesh::EquationSystems> XYCrossFieldSolver::_equation_systems

private

Holds the solve; kept past solve() because shape function evaluation needs it.

Definition at line 132 of file XYCrossFieldSolver.h.

Referenced by interpolatedCrossField(), solve(), and XYCrossFieldSolver().

_fe_type

const libMesh::FEType XYCrossFieldSolver::_fe_type

private

Linear Lagrange, the discretization of both components of z.

Definition at line 129 of file XYCrossFieldSolver.h.

Referenced by assembleLaplace(), interpolatedCrossField(), and XYCrossFieldSolver().

_imaginary_variable

unsigned int XYCrossFieldSolver::_imaginary_variable

private

Variable number of the imaginary part of z.

Definition at line 141 of file XYCrossFieldSolver.h.

Referenced by assembleLaplace(), extractNodalCrossField(), and XYCrossFieldSolver().

_mesh

libMesh::ReplicatedMesh XYCrossFieldSolver::_mesh

private

Private copy of the background triangulation, sharing the caller's node ids.

Definition at line 123 of file XYCrossFieldSolver.h.

Referenced by assembleLaplace(), extractNodalCrossField(), and XYCrossFieldSolver().

_nodal_cross_field

NodalCrossField XYCrossFieldSolver::_nodal_cross_field

private

Unit-magnitude nodal values of z, keyed by node id.

Definition at line 150 of file XYCrossFieldSolver.h.

Referenced by extractNodalCrossField(), interpolatedCrossField(), and nodalCrossField().

_point_locator

std::unique_ptr<libMesh::PointLocatorBase> XYCrossFieldSolver::_point_locator

private

Locates the element containing a query point.

Definition at line 144 of file XYCrossFieldSolver.h.

Referenced by XYCrossFieldSolver().

_real_variable

unsigned int XYCrossFieldSolver::_real_variable

private

Variable number of the real part of z.

Definition at line 138 of file XYCrossFieldSolver.h.

Referenced by assembleLaplace(), extractNodalCrossField(), and XYCrossFieldSolver().

_singular_nodes

std::vector<dof_id_type> XYCrossFieldSolver::_singular_nodes

private

Ascending ids of the nodes where z vanishes.

Definition at line 153 of file XYCrossFieldSolver.h.

Referenced by extractNodalCrossField(), and singularNodes().

_solved

bool XYCrossFieldSolver::_solved

private

Whether solve() has produced the results the queries read.

Definition at line 147 of file XYCrossFieldSolver.h.

Referenced by interpolatedCrossField(), and solve().

_system

libMesh::LinearImplicitSystem* XYCrossFieldSolver::_system

private

System holding both components of z; owned by _equation_systems.

Definition at line 135 of file XYCrossFieldSolver.h.

Referenced by assembleLaplace(), extractNodalCrossField(), solve(), and XYCrossFieldSolver().

---

The documentation for this class was generated from the following files:

• include/utils/XYCrossFieldSolver.h
• src/utils/XYCrossFieldSolver.C