<commit_before>// The libMesh Finite Element Library.
// Copyright (C) 2002-2012 <NAME>, <NAME>, <NAME>

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

// C++ Includes -------------------------------------
#include <set>
#include <algorithm> // for std::count, std::fill

// Local Includes -----------------------------------
#include "boundary_info.h" // needed for dirichlet constraints
#include "dense_matrix.h"
#include "dense_vector.h"
#include "dirichlet_boundaries.h"
#include "dof_map.h"
#include "elem.h"
#include "elem_range.h"
#include "fe_base.h"
#include "fe_interface.h"
#include "fe_type.h"
#include "libmesh_logging.h"
#include "system.h" // needed by enforce_constraints_exactly()
#include "mesh_base.h"
#include "numeric_vector.h" // for enforce_constraints_exactly()
#include "quadrature.h" // for dirichlet constraints
#include "parallel.h"
#include "periodic_boundaries.h"
#include "point_locator_base.h"
#include "threads.h"



// Anonymous namespace to hold helper classes
namespace {

using namespace libMesh;

  class ComputeConstraints
  {
  public:
    ComputeConstraints (DofConstraints &constraints,
			DofMap &dof_map,
#ifdef LIBMESH_ENABLE_PERIODIC
			PeriodicBoundaries &periodic_boundaries,
#endif
			const MeshBase &mesh,
			const unsigned int variable_number) :
      _constraints(constraints),
      _dof_map(dof_map),
#ifdef LIBMESH_ENABLE_PERIODIC
      _periodic_boundaries(periodic_boundaries),
#endif
      _mesh(mesh),
      _variable_number(variable_number)
    {}

    void operator()(const ConstElemRange &range) const
    {
      const Variable &var_description = _dof_map.variable(_variable_number);

#ifdef LIBMESH_ENABLE_PERIODIC
      AutoPtr<PointLocatorBase> point_locator;
      bool have_periodic_boundaries = !_periodic_boundaries.empty();
      if (have_periodic_boundaries)
        point_locator = _mesh.sub_point_locator();
#endif

      for (ConstElemRange::const_iterator it = range.begin(); it!=range.end(); ++it)
	if (var_description.active_on_subdomain((*it)->subdomain_id()))
	  {
#ifdef LIBMESH_ENABLE_AMR
	    FEInterface::compute_constraints (_constraints,
					      _dof_map,
					      _variable_number,
					      *it);
#endif
#ifdef LIBMESH_ENABLE_PERIODIC
	    // FIXME: periodic constraints won't work on a non-serial
	    // mesh unless it's kept ghost elements from opposing
	    // boundaries!
            if (have_periodic_boundaries)
	      FEInterface::compute_periodic_constraints (_constraints,
						         _dof_map,
						         _periodic_boundaries,
						         _mesh,
						         point_locator.get(),
						         _variable_number,
						         *it);
#endif
	  }
    }

  private:
    DofConstraints &_constraints;
    DofMap &_dof_map;
#ifdef LIBMESH_ENABLE_PERIODIC
    PeriodicBoundaries &_periodic_boundaries;
#endif
    const MeshBase &_mesh;
    const unsigned int _variable_number;
  };

#ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
  class ComputeNodeConstraints
  {
  public:
    ComputeNodeConstraints (NodeConstraints &node_constraints,
			    DofMap &dof_map,
#ifdef LIBMESH_ENABLE_PERIODIC
			    PeriodicBoundaries &periodic_boundaries,
#endif
			    const MeshBase &mesh) :
      _node_constraints(node_constraints),
      _dof_map(dof_map),
#ifdef LIBMESH_ENABLE_PERIODIC
      _periodic_boundaries(periodic_boundaries),
#endif
      _mesh(mesh)
    {}

    void operator()(const ConstElemRange &range) const
    {
#ifdef LIBMESH_ENABLE_PERIODIC
      AutoPtr<PointLocatorBase> point_locator;
      bool have_periodic_boundaries = !_periodic_boundaries.empty();
      if (have_periodic_boundaries)
        point_locator = _mesh.sub_point_locator();
#endif

      for (ConstElemRange::const_iterator it = range.begin(); it!=range.end(); ++it)
	{
#ifdef LIBMESH_ENABLE_AMR
	  FEBase::compute_node_constraints (_node_constraints, *it);
#endif
#ifdef LIBMESH_ENABLE_PERIODIC
	  // FIXME: periodic constraints won't work on a non-serial
	  // mesh unless it's kept ghost elements from opposing
	  // boundaries!
          if (have_periodic_boundaries)
	    FEBase::compute_periodic_node_constraints (_node_constraints,
						       _periodic_boundaries,
						       _mesh,
						       point_locator.get(),
						       *it);
#endif
	  }
    }

  private:
    NodeConstraints &_node_constraints;
    DofMap &_dof_map;
#ifdef LIBMESH_ENABLE_PERIODIC
    PeriodicBoundaries &_periodic_boundaries;
#endif
    const MeshBase &_mesh;
  };
#endif // LIBMESH_ENABLE_NODE_CONSTRAINTS


#ifdef LIBMESH_ENABLE_DIRICHLET

  /**
   * This class implements turning an arbitrary
   * boundary function into Dirichlet constraints.  It
   * may be executed in parallel on multiple threads.
   */
  class ConstrainDirichlet
  {
  private:
    DofMap                  &dof_map;
    const MeshBase          &mesh;
    const Real               time;
    const DirichletBoundary  dirichlet;

  public:
    ConstrainDirichlet (DofMap &dof_map_in,
		        const MeshBase &mesh_in,
		        const Real time_in,
		        const DirichletBoundary &dirichlet_in) :
    dof_map(dof_map_in),
    mesh(mesh_in),
    time(time_in),
    dirichlet(dirichlet_in) { }

    ConstrainDirichlet (const ConstrainDirichlet &in) :
    dof_map(in.dof_map),
    mesh(in.mesh),
    time(in.time),
    dirichlet(in.dirichlet) { }

    void operator()(const ConstElemRange &range) const
    {
      FunctionBase<Number> *f = dirichlet.f.get();
      FunctionBase<Gradient> *g = dirichlet.g.get();
      const std::set<boundary_id_type> &b = dirichlet.b;

      // We need data to project
      libmesh_assert(f);

      /**
       * This method examines an arbitrary boundary solution to calculate
       * corresponding Dirichlet constraints on the current mesh.  The
       * input function \p f gives the arbitrary solution.
       */

      // The dimensionality of the current mesh
      const unsigned int dim = mesh.mesh_dimension();

      // Boundary info for the current mesh
      const BoundaryInfo& boundary_info = *mesh.boundary_info;

      // The element matrix and RHS for projections.
      // Note that Ke is always real-valued, whereas
      // Fe may be complex valued if complex number
      // support is enabled
      DenseMatrix<Real> Ke;
      DenseVector<Number> Fe;
      // The new element coefficients
      DenseVector<Number> Ue;


      // Loop over all the variables we've been requested to project
      for (unsigned int v=0; v!=dirichlet.variables.size(); v++)
	{
	  const unsigned int var = dirichlet.variables[v];

	  const Variable& variable = dof_map.variable(var);

	  const FEType& fe_type = variable.type();

	  if (fe_type.family == SCALAR)
	    continue;

	  const unsigned int var_component =
	    variable.first_scalar_number();

	  // Get FE objects of the appropriate type
	  AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));

	  // Prepare variables for projection
	  AutoPtr<QBase> qedgerule (fe_type.default_quadrature_rule(1));
	  AutoPtr<QBase> qsiderule (fe_type.default_quadrature_rule(dim-1));

	  // The values of the shape functions at the quadrature
	  // points
	  const std::vector<std::vector<Real> >& phi = fe->get_phi();

	  // The gradients of the shape functions at the quadrature
	  // points on the child element.
	  const std::vector<std::vector<RealGradient> > *dphi = NULL;

	  const FEContinuity cont = fe->get_continuity();

	  if (cont == C_ONE)
	    {
	      // We'll need gradient data for a C1 projection
	      libmesh_assert(g);

	      const std::vector<std::vector<RealGradient> >&
		ref_dphi = fe->get_dphi();
	      dphi = &ref_dphi;
	    }

	  // The Jacobian * quadrature weight at the quadrature points
	  const std::vector<Real>& JxW =
	    fe->get_JxW();

	  // The XYZ locations of the quadrature points
	  const std::vector<Point>& xyz_values =
	    fe->get_xyz();

	  // The global DOF indices
	  std::vector<unsigned int> dof_indices;
	  // Side/edge DOF indices
	  std::vector<unsigned int> side_dofs;

	  // Iterate over all the elements in the range
	  for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it != range.end(); ++elem_it)
	    {
	      const Elem* elem = *elem_it;

	      // Per-subdomain variables don't need to be projected on
	      // elements where they're not active
	      if (!variable.active_on_subdomain(elem->subdomain_id()))
		continue;

	      // Find out which nodes, edges and sides are on a requested
	      // boundary:
	      std::vector<bool> is_boundary_node(elem->n_nodes(), false),
		is_boundary_edge(elem->n_edges(), false),
		is_boundary_side(elem->n_sides(), false);
	      for (unsigned char s=0; s != elem->n_sides(); ++s)
		{
		  // First see if this side has been requested
		  const std::vector<boundary_id_type>& bc_ids =
		    boundary_info.boundary_ids (elem, s);
		  bool do_this_side = false;
		  for (unsigned int i=0; i != bc_ids.size(); ++i)
		    if (b.count(bc_ids[i]))
		      {
			do_this_side = true;
			break;
		      }
		  if (!do_this_side)
		    continue;

		  // Then see what nodes and what edges are on it
		  for (unsigned int n=0; n != elem->n_nodes(); ++n)
		    if (elem->is_node_on_side(n,s))
		      is_boundary_node[n] = true;
		  for (unsigned int e=0; e != elem->n_edges(); ++e)
		    if (elem->is_edge_on_side(e,s))
		      is_boundary_edge[e] = true;
		}

	      // Update the DOF indices for this element based on
	      // the current mesh
	      dof_map.dof_indices (elem, dof_indices, var);

	      // The number of DOFs on the element
	      const unsigned int n_dofs = dof_indices.size();

	      // Fixed vs. free DoFs on edge/face projections
	      std::vector<char> dof_is_fixed(n_dofs, false); // bools
	      std::vector<int> free_dof(n_dofs, 0);

	      // The element type
	      const ElemType elem_type = elem->type();

	      // The number of nodes on the new element
	      const unsigned int n_nodes = elem->n_nodes();

	      // Zero the interpolated values
	      Ue.resize (n_dofs); Ue.zero();

	      // In general, we need a series of
	      // projections to ensure a unique and continuous
	      // solution.  We start by interpolating boundary nodes, then
	      // hold those fixed and project boundary edges, then hold
	      // those fixed and project boundary faces,

	      // Interpolate node values first
	      unsigned int current_dof = 0;
	      for (unsigned int n=0; n!= n_nodes; ++n)
		{
		  // FIXME: this should go through the DofMap,
		  // not duplicate dof_indices code badly!
		  const unsigned int nc =
		    FEInterface::n_dofs_at_node (dim, fe_type, elem_type,
						 n);
		  if (!elem->is_vertex(n) || !is_boundary_node[n])
		    {
		      current_dof += nc;
		      continue;
		    }
		  if (cont == DISCONTINUOUS)
		    {
		      libmesh_assert(nc == 0);
		    }
		  // Assume that C_ZERO elements have a single nodal
		  // value shape function
		  else if (cont == C_ZERO)
		    {
		      libmesh_assert(nc == 1);
		      Ue(current_dof) =
                        f->component(var_component,
                                     elem->point(n), time);
		      dof_is_fixed[current_dof] = true;
		      current_dof++;
		    }
		  // The hermite element vertex shape functions are weird
		  else if (fe_type.family == HERMITE)
		    {
		      Ue(current_dof) =
                        f->component(var_component,
                                     elem->point(n), time);
		      dof_is_fixed[current_dof] = true;
		      current_dof++;
		      Gradient grad =
                        g->component(var_component,
                                     elem->point(n), time);
		      // x derivative
		      Ue(current_dof) = grad(0);
		      dof_is_fixed[current_dof] = true;
		      current_dof++;
		      if (dim > 1)
			{
			  // We'll finite difference mixed derivatives
			  Point nxminus = elem->point(n),
                            nxplus = elem->point(n);
			  nxminus(0) -= TOLERANCE;
			  nxplus(0) += TOLERANCE;
			  Gradient gxminus =
                            g->component(var_component,
                                         nxminus, time);
			  Gradient gxplus =
                            g->component(var_component,
                                         nxplus, time);
			  // y derivative
			  Ue(current_dof) = grad(1);
			  dof_is_fixed[current_dof] = true;
			  current_dof++;
			  // xy derivative
			  Ue(current_dof) = (gxplus(1) - gxminus(1))
			    / 2. / TOLERANCE;
			  dof_is_fixed[current_dof] = true;
			  current_dof++;

			  if (dim > 2)
			    {
			      // z derivative
			      Ue(current_dof) = grad(2);
			      dof_is_fixed[current_dof] = true;
			      current_dof++;
			      // xz derivative
			      Ue(current_dof) = (gxplus(2) - gxminus(2))
				/ 2. / TOLERANCE;
			      dof_is_fixed[current_dof] = true;
			      current_dof++;
			      // We need new points for yz
			      Point nyminus = elem->point(n),
                                nyplus = elem->point(n);
			      nyminus(1) -= TOLERANCE;
			      nyplus(1) += TOLERANCE;
			      Gradient gyminus =
                                g->component(var_component,
                                             nyminus, time);
			      Gradient gyplus =
                                g->component(var_component,
                                             nyplus, time);
			      // xz derivative
			      Ue(current_dof) = (gyplus(2) - gyminus(2))
				/ 2. / TOLERANCE;
			      dof_is_fixed[current_dof] = true;
			      current_dof++;
			      // Getting a 2nd order xyz is more tedious
			      Point nxmym = elem->point(n),
                                nxmyp = elem->point(n),
                                nxpym = elem->point(n),
                                nxpyp = elem->point(n);
			      nxmym(0) -= TOLERANCE;
			      nxmym(1) -= TOLERANCE;
			      nxmyp(0) -= TOLERANCE;
			      nxmyp(1) += TOLERANCE;
			      nxpym(0) += TOLERANCE;
			      nxpym(1) -= TOLERANCE;
			      nxpyp(0) += TOLERANCE;
			      nxpyp(1) += TOLERANCE;
			      Gradient gxmym =
                                g->component(var_component,
                                             nxmym, time);
			      Gradient gxmyp =
                                g->component(var_component,
                                             nxmyp, time);
			      Gradient gxpym =
                                g->component(var_component,
                                             nxpym, time);
			      Gradient gxpyp =
                                g->component(var_component,
                                             nxpyp, time);
			      Number gxzplus = (gxpyp(2) - gxmyp(2))
				/ 2. / TOLERANCE;
			      Number gxzminus = (gxpym(2) - gxmym(2))
				/ 2. / TOLERANCE;
			      // xyz derivative
			      Ue(current_dof) = (gxzplus - gxzminus)
				/ 2. / TOLERANCE;
			      dof_is_fixed[current_dof] = true;
			      current_dof++;
			    }
			}
		    }
		  // Assume that other C_ONE elements have a single nodal
		  // value shape function and nodal gradient component
		  // shape functions
		  else if (cont == C_ONE)
		    {
		      libmesh_assert(nc == 1 + dim);
		      Ue(current_dof) =
                        f->component(var_component,
                                     elem->point(n), time);
		      dof_is_fixed[current_dof] = true;
		      current_dof++;
		      Gradient grad =
                        g->component(var_component,
                                     elem->point(n), time);
		      for (unsigned int i=0; i!= dim; ++i)
			{
			  Ue(current_dof) = grad(i);
			  dof_is_fixed[current_dof] = true;
			  current_dof++;
			}
		    }
		  else
		    libmesh_error();
		}

	      // In 3D, project any edge values next
	      if (dim > 2 && cont != DISCONTINUOUS)
		for (unsigned int e=0; e != elem->n_edges(); ++e)
		  {
		    if (!is_boundary_edge[e])
		      continue;

		    FEInterface::dofs_on_edge(elem, dim, fe_type, e,
					      side_dofs);

		    // Some edge dofs are on nodes and already
		    // fixed, others are free to calculate
		    unsigned int free_dofs = 0;
		    for (unsigned int i=0; i != side_dofs.size(); ++i)
		      if (!dof_is_fixed[side_dofs[i]])
			free_dof[free_dofs++] = i;

		    // There may be nothing to project
		    if (!free_dofs)
		      continue;

		    Ke.resize (free_dofs, free_dofs); Ke.zero();
		    Fe.resize (free_dofs); Fe.zero();
		    // The new edge coefficients
		    DenseVector<Number> Uedge(free_dofs);

		    // Initialize FE data on the edge
		    fe->attach_quadrature_rule (qedgerule.get());
		    fe->edge_reinit (elem, e);
		    const unsigned int n_qp = qedgerule->n_points();

		    // Loop over the quadrature points
		    for (unsigned int qp=0; qp<n_qp; qp++)
		      {
			// solution at the quadrature point
			Number fineval =
                          f->component(var_component,
                                       xyz_values[qp], time);
			// solution grad at the quadrature point
			Gradient finegrad;
			if (cont == C_ONE)
			  finegrad =
                            g->component(var_component,
                                         xyz_values[qp], time);

			// Form edge projection matrix
			for (unsigned int sidei=0, freei=0;
			     sidei != side_dofs.size(); ++sidei)
			  {
			    unsigned int i = side_dofs[sidei];
			    // fixed DoFs aren't test functions
			    if (dof_is_fixed[i])
			      continue;
			    for (unsigned int sidej=0, freej=0;
				 sidej != side_dofs.size(); ++sidej)
			      {
				unsigned int j = side_dofs[sidej];
				if (dof_is_fixed[j])
				  Fe(freei) -= phi[i][qp] * phi[j][qp] *
				    JxW[qp] * Ue(j);
				else
				  Ke(freei,freej) += phi[i][qp] *
				    phi[j][qp] * JxW[qp];
				if (cont == C_ONE)
				  {
				    if (dof_is_fixed[j])
				      Fe(freei) -= ((*dphi)[i][qp] *
						    (*dphi)[j][qp]) *
					JxW[qp] * Ue(j);
				    else
				      Ke(freei,freej) += ((*dphi)[i][qp] *
							  (*dphi)[j][qp])
					* JxW[qp];
				  }
				if (!dof_is_fixed[j])
				  freej++;
			      }
			    Fe(freei) += phi[i][qp] * fineval * JxW[qp];
			    if (cont == C_ONE)
			      Fe(freei) += (finegrad * (*dphi)[i][qp]) *
				JxW[qp];
			    freei++;
			  }
		      }

		    Ke.cholesky_solve(Fe, Uedge);

		    // Transfer new edge solutions to element
		    for (unsigned int i=0; i != free_dofs; ++i)
		      {
			Number &ui = Ue(side_dofs[free_dof[i]]);
			libmesh_assert(std::abs(ui) < TOLERANCE ||
				       std::abs(ui - Uedge(i)) < TOLERANCE);
			ui = Uedge(i);
			dof_is_fixed[side_dofs[free_dof[i]]] = true;
		      }
		  }

	      // Project any side values (edges in 2D, faces in 3D)
	      if (dim > 1 && cont != DISCONTINUOUS)
		for (unsigned int s=0; s != elem->n_sides(); ++s)
		  {
		    if (!is_boundary_side[s])
		      continue;

		    FEInterface::dofs_on_side(elem, dim, fe_type, s,
					      side_dofs);

		    // Some side dofs are on nodes/edges and already
		    // fixed, others are free to calculate
		    unsigned int free_dofs = 0;
		    for (unsigned int i=0; i != side_dofs.size(); ++i)
		      if (!dof_is_fixed[side_dofs[i]])
			free_dof[free_dofs++] = i;

		    // There may be nothing to project
		    if (!free_dofs)
		      continue;

		    Ke.resize (free_dofs, free_dofs); Ke.zero();
		    Fe.resize (free_dofs); Fe.zero();
		    // The new side coefficients
		    DenseVector<Number> Uside(free_dofs);

		    // Initialize FE data on the side
		    fe->attach_quadrature_rule (qsiderule.get());
		    fe->reinit (elem, s);
		    const unsigned int n_qp = qsiderule->n_points();

		    // Loop over the quadrature points
		    for (unsigned int qp=0; qp<n_qp; qp++)
		      {
			// solution at the quadrature point
			Number fineval =
                          f->component(var_component,
                                       xyz_values[qp], time);
			// solution grad at the quadrature point
			Gradient finegrad;
			if (cont == C_ONE)
			  finegrad =
                            g->component(var_component,
                                         xyz_values[qp], time);

			// Form side projection matrix
			for (unsigned int sidei=0, freei=0;
			     sidei != side_dofs.size(); ++sidei)
			  {
			    unsigned int i = side_dofs[sidei];
			    // fixed DoFs aren't test functions
			    if (dof_is_fixed[i])
			      continue;
			    for (unsigned int sidej=0, freej=0;
				 sidej != side_dofs.size(); ++sidej)
			      {
				unsigned int j = side_dofs[sidej];
				if (dof_is_fixed[j])
				  Fe(freei) -= phi[i][qp] * phi[j][qp] *
				    JxW[qp] * Ue(j);
				else
				  Ke(freei,freej) += phi[i][qp] *
				    phi[j][qp] * JxW[qp];
				if (cont == C_ONE)
				  {
				    if (dof_is_fixed[j])
				      Fe(freei) -= ((*dphi)[i][qp] *
						    (*dphi)[j][qp]) *
					JxW[qp] * Ue(j);
				    else
				      Ke(freei,freej) += ((*dphi)[i][qp] *
							  (*dphi)[j][qp])
					* JxW[qp];
				  }
				if (!dof_is_fixed[j])
				  freej++;
			      }
			    Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
			    if (cont == C_ONE)
			      Fe(freei) += (finegrad * (*dphi)[i][qp]) *
				JxW[qp];
			    freei++;
			  }
		      }

		    Ke.cholesky_solve(Fe, Uside);

		    // Transfer new side solutions to element
		    for (unsigned int i=0; i != free_dofs; ++i)
		      {
			Number &ui = Ue(side_dofs[free_dof[i]]);
			libmesh_assert(std::abs(ui) < TOLERANCE ||
				       std::abs(ui - Uside(i)) < TOLERANCE);
			ui = Uside(i);
			dof_is_fixed[side_dofs[free_dof[i]]] = true;
		      }
		  }

	      // Lock the DofConstraints since it is shared among threads.
	      {
		Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);

		for (unsigned int i = 0; i < n_dofs; i++)
		  {
		    DofConstraintRow empty_row;
		    dof_map.add_constraint_row
		      (dof_indices[i], empty_row, 
		       Ue(i), /* forbid_constraint_overwrite = */ true);
		  }
	      }
	    }  // end elem loop
	} // end variables loop
    }
  };

#endif // LIBMESH_ENABLE_DIRICHLET


} // anonymous namespace



namespace libMesh
{

// ------------------------------------------------------------
// DofMap member functions

#if defined(LIBMESH_ENABLE_AMR) || \
    defined(LIBMESH_ENABLE_PERIODIC) || \
    defined(LIBMESH_ENABLE_DIRICHLET)

void DofMap::create_dof_constraints(const MeshBase& mesh, Real time)
{
  START_LOG("create_dof_constraints()", "DofMap");
<commit_msg>Don't try to add Dirichlet constraints on already constrained DoFs.<commit_after>// The libMesh Finite Element Library.
// Copyright (C) 2002-2012 <NAME>, <NAME>, <NAME>

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

// C++ Includes -------------------------------------
#include <set>
#include <algorithm> // for std::count, std::fill

// Local Includes -----------------------------------
#include "boundary_info.h" // needed for dirichlet constraints
#include "dense_matrix.h"
#include "dense_vector.h"
#include "dirichlet_boundaries.h"
#include "dof_map.h"
#include "elem.h"
#include "elem_range.h"
#include "fe_base.h"
#include "fe_interface.h"
#include "fe_type.h"
#include "libmesh_logging.h"
#include "system.h" // needed by enforce_constraints_exactly()
#include "mesh_base.h"
#include "numeric_vector.h" // for enforce_constraints_exactly()
#include "quadrature.h" // for dirichlet constraints
#include "parallel.h"
#include "periodic_boundaries.h"
#include "point_locator_base.h"
#include "threads.h"



// Anonymous namespace to hold helper classes
namespace {

using namespace libMesh;

  class ComputeConstraints
  {
  public:
    ComputeConstraints (DofConstraints &constraints,
			DofMap &dof_map,
#ifdef LIBMESH_ENABLE_PERIODIC
			PeriodicBoundaries &periodic_boundaries,
#endif
			const MeshBase &mesh,
			const unsigned int variable_number) :
      _constraints(constraints),
      _dof_map(dof_map),
#ifdef LIBMESH_ENABLE_PERIODIC
      _periodic_boundaries(periodic_boundaries),
#endif
      _mesh(mesh),
      _variable_number(variable_number)
    {}

    void operator()(const ConstElemRange &range) const
    {
      const Variable &var_description = _dof_map.variable(_variable_number);

#ifdef LIBMESH_ENABLE_PERIODIC
      AutoPtr<PointLocatorBase> point_locator;
      bool have_periodic_boundaries = !_periodic_boundaries.empty();
      if (have_periodic_boundaries)
        point_locator = _mesh.sub_point_locator();
#endif

      for (ConstElemRange::const_iterator it = range.begin(); it!=range.end(); ++it)
	if (var_description.active_on_subdomain((*it)->subdomain_id()))
	  {
#ifdef LIBMESH_ENABLE_AMR
	    FEInterface::compute_constraints (_constraints,
					      _dof_map,
					      _variable_number,
					      *it);
#endif
#ifdef LIBMESH_ENABLE_PERIODIC
	    // FIXME: periodic constraints won't work on a non-serial
	    // mesh unless it's kept ghost elements from opposing
	    // boundaries!
            if (have_periodic_boundaries)
	      FEInterface::compute_periodic_constraints (_constraints,
						         _dof_map,
						         _periodic_boundaries,
						         _mesh,
						         point_locator.get(),
						         _variable_number,
						         *it);
#endif
	  }
    }

  private:
    DofConstraints &_constraints;
    DofMap &_dof_map;
#ifdef LIBMESH_ENABLE_PERIODIC
    PeriodicBoundaries &_periodic_boundaries;
#endif
    const MeshBase &_mesh;
    const unsigned int _variable_number;
  };

#ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
  class ComputeNodeConstraints
  {
  public:
    ComputeNodeConstraints (NodeConstraints &node_constraints,
			    DofMap &dof_map,
#ifdef LIBMESH_ENABLE_PERIODIC
			    PeriodicBoundaries &periodic_boundaries,
#endif
			    const MeshBase &mesh) :
      _node_constraints(node_constraints),
      _dof_map(dof_map),
#ifdef LIBMESH_ENABLE_PERIODIC
      _periodic_boundaries(periodic_boundaries),
#endif
      _mesh(mesh)
    {}

    void operator()(const ConstElemRange &range) const
    {
#ifdef LIBMESH_ENABLE_PERIODIC
      AutoPtr<PointLocatorBase> point_locator;
      bool have_periodic_boundaries = !_periodic_boundaries.empty();
      if (have_periodic_boundaries)
        point_locator = _mesh.sub_point_locator();
#endif

      for (ConstElemRange::const_iterator it = range.begin(); it!=range.end(); ++it)
	{
#ifdef LIBMESH_ENABLE_AMR
	  FEBase::compute_node_constraints (_node_constraints, *it);
#endif
#ifdef LIBMESH_ENABLE_PERIODIC
	  // FIXME: periodic constraints won't work on a non-serial
	  // mesh unless it's kept ghost elements from opposing
	  // boundaries!
          if (have_periodic_boundaries)
	    FEBase::compute_periodic_node_constraints (_node_constraints,
						       _periodic_boundaries,
						       _mesh,
						       point_locator.get(),
						       *it);
#endif
	  }
    }

  private:
    NodeConstraints &_node_constraints;
    DofMap &_dof_map;
#ifdef LIBMESH_ENABLE_PERIODIC
    PeriodicBoundaries &_periodic_boundaries;
#endif
    const MeshBase &_mesh;
  };
#endif // LIBMESH_ENABLE_NODE_CONSTRAINTS


#ifdef LIBMESH_ENABLE_DIRICHLET

  /**
   * This class implements turning an arbitrary
   * boundary function into Dirichlet constraints.  It
   * may be executed in parallel on multiple threads.
   */
  class ConstrainDirichlet
  {
  private:
    DofMap                  &dof_map;
    const MeshBase          &mesh;
    const Real               time;
    const DirichletBoundary  dirichlet;

  public:
    ConstrainDirichlet (DofMap &dof_map_in,
		        const MeshBase &mesh_in,
		        const Real time_in,
		        const DirichletBoundary &dirichlet_in) :
    dof_map(dof_map_in),
    mesh(mesh_in),
    time(time_in),
    dirichlet(dirichlet_in) { }

    ConstrainDirichlet (const ConstrainDirichlet &in) :
    dof_map(in.dof_map),
    mesh(in.mesh),
    time(in.time),
    dirichlet(in.dirichlet) { }

    void operator()(const ConstElemRange &range) const
    {
      FunctionBase<Number> *f = dirichlet.f.get();
      FunctionBase<Gradient> *g = dirichlet.g.get();
      const std::set<boundary_id_type> &b = dirichlet.b;

      // We need data to project
      libmesh_assert(f);

      /**
       * This method examines an arbitrary boundary solution to calculate
       * corresponding Dirichlet constraints on the current mesh.  The
       * input function \p f gives the arbitrary solution.
       */

      // The dimensionality of the current mesh
      const unsigned int dim = mesh.mesh_dimension();

      // Boundary info for the current mesh
      const BoundaryInfo& boundary_info = *mesh.boundary_info;

      // The element matrix and RHS for projections.
      // Note that Ke is always real-valued, whereas
      // Fe may be complex valued if complex number
      // support is enabled
      DenseMatrix<Real> Ke;
      DenseVector<Number> Fe;
      // The new element coefficients
      DenseVector<Number> Ue;


      // Loop over all the variables we've been requested to project
      for (unsigned int v=0; v!=dirichlet.variables.size(); v++)
	{
	  const unsigned int var = dirichlet.variables[v];

	  const Variable& variable = dof_map.variable(var);

	  const FEType& fe_type = variable.type();

	  if (fe_type.family == SCALAR)
	    continue;

	  const unsigned int var_component =
	    variable.first_scalar_number();

	  // Get FE objects of the appropriate type
	  AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));

	  // Prepare variables for projection
	  AutoPtr<QBase> qedgerule (fe_type.default_quadrature_rule(1));
	  AutoPtr<QBase> qsiderule (fe_type.default_quadrature_rule(dim-1));

	  // The values of the shape functions at the quadrature
	  // points
	  const std::vector<std::vector<Real> >& phi = fe->get_phi();

	  // The gradients of the shape functions at the quadrature
	  // points on the child element.
	  const std::vector<std::vector<RealGradient> > *dphi = NULL;

	  const FEContinuity cont = fe->get_continuity();

	  if (cont == C_ONE)
	    {
	      // We'll need gradient data for a C1 projection
	      libmesh_assert(g);

	      const std::vector<std::vector<RealGradient> >&
		ref_dphi = fe->get_dphi();
	      dphi = &ref_dphi;
	    }

	  // The Jacobian * quadrature weight at the quadrature points
	  const std::vector<Real>& JxW =
	    fe->get_JxW();

	  // The XYZ locations of the quadrature points
	  const std::vector<Point>& xyz_values =
	    fe->get_xyz();

	  // The global DOF indices
	  std::vector<unsigned int> dof_indices;
	  // Side/edge DOF indices
	  std::vector<unsigned int> side_dofs;

	  // Iterate over all the elements in the range
	  for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it != range.end(); ++elem_it)
	    {
	      const Elem* elem = *elem_it;

	      // Per-subdomain variables don't need to be projected on
	      // elements where they're not active
	      if (!variable.active_on_subdomain(elem->subdomain_id()))
		continue;

	      // Find out which nodes, edges and sides are on a requested
	      // boundary:
	      std::vector<bool> is_boundary_node(elem->n_nodes(), false),
		is_boundary_edge(elem->n_edges(), false),
		is_boundary_side(elem->n_sides(), false);
	      for (unsigned char s=0; s != elem->n_sides(); ++s)
		{
		  // First see if this side has been requested
		  const std::vector<boundary_id_type>& bc_ids =
		    boundary_info.boundary_ids (elem, s);
		  bool do_this_side = false;
		  for (unsigned int i=0; i != bc_ids.size(); ++i)
		    if (b.count(bc_ids[i]))
		      {
			do_this_side = true;
			break;
		      }
		  if (!do_this_side)
		    continue;

		  // Then see what nodes and what edges are on it
		  for (unsigned int n=0; n != elem->n_nodes(); ++n)
		    if (elem->is_node_on_side(n,s))
		      is_boundary_node[n] = true;
		  for (unsigned int e=0; e != elem->n_edges(); ++e)
		    if (elem->is_edge_on_side(e,s))
		      is_boundary_edge[e] = true;
		}

	      // Update the DOF indices for this element based on
	      // the current mesh
	      dof_map.dof_indices (elem, dof_indices, var);

	      // The number of DOFs on the element
	      const unsigned int n_dofs = dof_indices.size();

	      // Fixed vs. free DoFs on edge/face projections
	      std::vector<char> dof_is_fixed(n_dofs, false); // bools
	      std::vector<int> free_dof(n_dofs, 0);

	      // The element type
	      const ElemType elem_type = elem->type();

	      // The number of nodes on the new element
	      const unsigned int n_nodes = elem->n_nodes();

	      // Zero the interpolated values
	      Ue.resize (n_dofs); Ue.zero();

	      // In general, we need a series of
	      // projections to ensure a unique and continuous
	      // solution.  We start by interpolating boundary nodes, then
	      // hold those fixed and project boundary edges, then hold
	      // those fixed and project boundary faces,

	      // Interpolate node values first
	      unsigned int current_dof = 0;
	      for (unsigned int n=0; n!= n_nodes; ++n)
		{
		  // FIXME: this should go through the DofMap,
		  // not duplicate dof_indices code badly!
		  const unsigned int nc =
		    FEInterface::n_dofs_at_node (dim, fe_type, elem_type,
						 n);
		  if (!elem->is_vertex(n) || !is_boundary_node[n])
		    {
		      current_dof += nc;
		      continue;
		    }
		  if (cont == DISCONTINUOUS)
		    {
		      libmesh_assert(nc == 0);
		    }
		  // Assume that C_ZERO elements have a single nodal
		  // value shape function
		  else if (cont == C_ZERO)
		    {
		      libmesh_assert(nc == 1);
		      Ue(current_dof) =
                        f->component(var_component,
                                     elem->point(n), time);
		      dof_is_fixed[current_dof] = true;
		      current_dof++;
		    }
		  // The hermite element vertex shape functions are weird
		  else if (fe_type.family == HERMITE)
		    {
		      Ue(current_dof) =
                        f->component(var_component,
                                     elem->point(n), time);
		      dof_is_fixed[current_dof] = true;
		      current_dof++;
		      Gradient grad =
                        g->component(var_component,
                                     elem->point(n), time);
		      // x derivative
		      Ue(current_dof) = grad(0);
		      dof_is_fixed[current_dof] = true;
		      current_dof++;
		      if (dim > 1)
			{
			  // We'll finite difference mixed derivatives
			  Point nxminus = elem->point(n),
                            nxplus = elem->point(n);
			  nxminus(0) -= TOLERANCE;
			  nxplus(0) += TOLERANCE;
			  Gradient gxminus =
                            g->component(var_component,
                                         nxminus, time);
			  Gradient gxplus =
                            g->component(var_component,
                                         nxplus, time);
			  // y derivative
			  Ue(current_dof) = grad(1);
			  dof_is_fixed[current_dof] = true;
			  current_dof++;
			  // xy derivative
			  Ue(current_dof) = (gxplus(1) - gxminus(1))
			    / 2. / TOLERANCE;
			  dof_is_fixed[current_dof] = true;
			  current_dof++;

			  if (dim > 2)
			    {
			      // z derivative
			      Ue(current_dof) = grad(2);
			      dof_is_fixed[current_dof] = true;
			      current_dof++;
			      // xz derivative
			      Ue(current_dof) = (gxplus(2) - gxminus(2))
				/ 2. / TOLERANCE;
			      dof_is_fixed[current_dof] = true;
			      current_dof++;
			      // We need new points for yz
			      Point nyminus = elem->point(n),
                                nyplus = elem->point(n);
			      nyminus(1) -= TOLERANCE;
			      nyplus(1) += TOLERANCE;
			      Gradient gyminus =
                                g->component(var_component,
                                             nyminus, time);
			      Gradient gyplus =
                                g->component(var_component,
                                             nyplus, time);
			      // xz derivative
			      Ue(current_dof) = (gyplus(2) - gyminus(2))
				/ 2. / TOLERANCE;
			      dof_is_fixed[current_dof] = true;
			      current_dof++;
			      // Getting a 2nd order xyz is more tedious
			      Point nxmym = elem->point(n),
                                nxmyp = elem->point(n),
                                nxpym = elem->point(n),
                                nxpyp = elem->point(n);
			      nxmym(0) -= TOLERANCE;
			      nxmym(1) -= TOLERANCE;
			      nxmyp(0) -= TOLERANCE;
			      nxmyp(1) += TOLERANCE;
			      nxpym(0) += TOLERANCE;
			      nxpym(1) -= TOLERANCE;
			      nxpyp(0) += TOLERANCE;
			      nxpyp(1) += TOLERANCE;
			      Gradient gxmym =
                                g->component(var_component,
                                             nxmym, time);
			      Gradient gxmyp =
                                g->component(var_component,
                                             nxmyp, time);
			      Gradient gxpym =
                                g->component(var_component,
                                             nxpym, time);
			      Gradient gxpyp =
                                g->component(var_component,
                                             nxpyp, time);
			      Number gxzplus = (gxpyp(2) - gxmyp(2))
				/ 2. / TOLERANCE;
			      Number gxzminus = (gxpym(2) - gxmym(2))
				/ 2. / TOLERANCE;
			      // xyz derivative
			      Ue(current_dof) = (gxzplus - gxzminus)
				/ 2. / TOLERANCE;
			      dof_is_fixed[current_dof] = true;
			      current_dof++;
			    }
			}
		    }
		  // Assume that other C_ONE elements have a single nodal
		  // value shape function and nodal gradient component
		  // shape functions
		  else if (cont == C_ONE)
		    {
		      libmesh_assert(nc == 1 + dim);
		      Ue(current_dof) =
                        f->component(var_component,
                                     elem->point(n), time);
		      dof_is_fixed[current_dof] = true;
		      current_dof++;
		      Gradient grad =
                        g->component(var_component,
                                     elem->point(n), time);
		      for (unsigned int i=0; i!= dim; ++i)
			{
			  Ue(current_dof) = grad(i);
			  dof_is_fixed[current_dof] = true;
			  current_dof++;
			}
		    }
		  else
		    libmesh_error();
		}

	      // In 3D, project any edge values next
	      if (dim > 2 && cont != DISCONTINUOUS)
		for (unsigned int e=0; e != elem->n_edges(); ++e)
		  {
		    if (!is_boundary_edge[e])
		      continue;

		    FEInterface::dofs_on_edge(elem, dim, fe_type, e,
					      side_dofs);

		    // Some edge dofs are on nodes and already
		    // fixed, others are free to calculate
		    unsigned int free_dofs = 0;
		    for (unsigned int i=0; i != side_dofs.size(); ++i)
		      if (!dof_is_fixed[side_dofs[i]])
			free_dof[free_dofs++] = i;

		    // There may be nothing to project
		    if (!free_dofs)
		      continue;

		    Ke.resize (free_dofs, free_dofs); Ke.zero();
		    Fe.resize (free_dofs); Fe.zero();
		    // The new edge coefficients
		    DenseVector<Number> Uedge(free_dofs);

		    // Initialize FE data on the edge
		    fe->attach_quadrature_rule (qedgerule.get());
		    fe->edge_reinit (elem, e);
		    const unsigned int n_qp = qedgerule->n_points();

		    // Loop over the quadrature points
		    for (unsigned int qp=0; qp<n_qp; qp++)
		      {
			// solution at the quadrature point
			Number fineval =
                          f->component(var_component,
                                       xyz_values[qp], time);
			// solution grad at the quadrature point
			Gradient finegrad;
			if (cont == C_ONE)
			  finegrad =
                            g->component(var_component,
                                         xyz_values[qp], time);

			// Form edge projection matrix
			for (unsigned int sidei=0, freei=0;
			     sidei != side_dofs.size(); ++sidei)
			  {
			    unsigned int i = side_dofs[sidei];
			    // fixed DoFs aren't test functions
			    if (dof_is_fixed[i])
			      continue;
			    for (unsigned int sidej=0, freej=0;
				 sidej != side_dofs.size(); ++sidej)
			      {
				unsigned int j = side_dofs[sidej];
				if (dof_is_fixed[j])
				  Fe(freei) -= phi[i][qp] * phi[j][qp] *
				    JxW[qp] * Ue(j);
				else
				  Ke(freei,freej) += phi[i][qp] *
				    phi[j][qp] * JxW[qp];
				if (cont == C_ONE)
				  {
				    if (dof_is_fixed[j])
				      Fe(freei) -= ((*dphi)[i][qp] *
						    (*dphi)[j][qp]) *
					JxW[qp] * Ue(j);
				    else
				      Ke(freei,freej) += ((*dphi)[i][qp] *
							  (*dphi)[j][qp])
					* JxW[qp];
				  }
				if (!dof_is_fixed[j])
				  freej++;
			      }
			    Fe(freei) += phi[i][qp] * fineval * JxW[qp];
			    if (cont == C_ONE)
			      Fe(freei) += (finegrad * (*dphi)[i][qp]) *
				JxW[qp];
			    freei++;
			  }
		      }

		    Ke.cholesky_solve(Fe, Uedge);

		    // Transfer new edge solutions to element
		    for (unsigned int i=0; i != free_dofs; ++i)
		      {
			Number &ui = Ue(side_dofs[free_dof[i]]);
			libmesh_assert(std::abs(ui) < TOLERANCE ||
				       std::abs(ui - Uedge(i)) < TOLERANCE);
			ui = Uedge(i);
			dof_is_fixed[side_dofs[free_dof[i]]] = true;
		      }
		  }

	      // Project any side values (edges in 2D, faces in 3D)
	      if (dim > 1 && cont != DISCONTINUOUS)
		for (unsigned int s=0; s != elem->n_sides(); ++s)
		  {
		    if (!is_boundary_side[s])
		      continue;

		    FEInterface::dofs_on_side(elem, dim, fe_type, s,
					      side_dofs);

		    // Some side dofs are on nodes/edges and already
		    // fixed, others are free to calculate
		    unsigned int free_dofs = 0;
		    for (unsigned int i=0; i != side_dofs.size(); ++i)
		      if (!dof_is_fixed[side_dofs[i]])
			free_dof[free_dofs++] = i;

		    // There may be nothing to project
		    if (!free_dofs)
		      continue;

		    Ke.resize (free_dofs, free_dofs); Ke.zero();
		    Fe.resize (free_dofs); Fe.zero();
		    // The new side coefficients
		    DenseVector<Number> Uside(free_dofs);

		    // Initialize FE data on the side
		    fe->attach_quadrature_rule (qsiderule.get());
		    fe->reinit (elem, s);
		    const unsigned int n_qp = qsiderule->n_points();

		    // Loop over the quadrature points
		    for (unsigned int qp=0; qp<n_qp; qp++)
		      {
			// solution at the quadrature point
			Number fineval =
                          f->component(var_component,
                                       xyz_values[qp], time);
			// solution grad at the quadrature point
			Gradient finegrad;
			if (cont == C_ONE)
			  finegrad =
                            g->component(var_component,
                                         xyz_values[qp], time);

			// Form side projection matrix
			for (unsigned int sidei=0, freei=0;
			     sidei != side_dofs.size(); ++sidei)
			  {
			    unsigned int i = side_dofs[sidei];
			    // fixed DoFs aren't test functions
			    if (dof_is_fixed[i])
			      continue;
			    for (unsigned int sidej=0, freej=0;
				 sidej != side_dofs.size(); ++sidej)
			      {
				unsigned int j = side_dofs[sidej];
				if (dof_is_fixed[j])
				  Fe(freei) -= phi[i][qp] * phi[j][qp] *
				    JxW[qp] * Ue(j);
				else
				  Ke(freei,freej) += phi[i][qp] *
				    phi[j][qp] * JxW[qp];
				if (cont == C_ONE)
				  {
				    if (dof_is_fixed[j])
				      Fe(freei) -= ((*dphi)[i][qp] *
						    (*dphi)[j][qp]) *
					JxW[qp] * Ue(j);
				    else
				      Ke(freei,freej) += ((*dphi)[i][qp] *
							  (*dphi)[j][qp])
					* JxW[qp];
				  }
				if (!dof_is_fixed[j])
				  freej++;
			      }
			    Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
			    if (cont == C_ONE)
			      Fe(freei) += (finegrad * (*dphi)[i][qp]) *
				JxW[qp];
			    freei++;
			  }
		      }

		    Ke.cholesky_solve(Fe, Uside);

		    // Transfer new side solutions to element
		    for (unsigned int i=0; i != free_dofs; ++i)
		      {
			Number &ui = Ue(side_dofs[free_dof[i]]);
			libmesh_assert(std::abs(ui) < TOLERANCE ||
				       std::abs(ui - Uside(i)) < TOLERANCE);
			ui = Uside(i);
			dof_is_fixed[side_dofs[free_dof[i]]] = true;
		      }
		  }

	      // Lock the DofConstraints since it is shared among threads.
	      {
		Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);

		for (unsigned int i = 0; i < n_dofs; i++)
		  {
		    DofConstraintRow empty_row;
		    if (!dof_map.is_constrained_dof(dof_indices[i]))
		      dof_map.add_constraint_row
		        (dof_indices[i], empty_row, 
		         Ue(i), /* forbid_constraint_overwrite = */ true);
		  }
	      }
	    }  // end elem loop
	} // end variables loop
    }
  };

#endif // LIBMESH_ENABLE_DIRICHLET


} // anonymous namespace



namespace libMesh
{

// ------------------------------------------------------------
// DofMap member functions

#if defined(LIBMESH_ENABLE_AMR) || \
    defined(LIBMESH_ENABLE_PERIODIC) || \
    defined(LIBMESH_ENABLE_DIRICHLET)

void DofMap::create_dof_constraints(const MeshBase& mesh, Real time)
{
  START_LOG("create_dof_constraints()", "DofMap");

{'fraction_non_alphanumeric': 0.11706914165235725, 'fraction_numerical': 0.004216749625057056, 'mean_word_length': 2.085714285714286, 'pattern_counts': {'":': 0, '<': 81, '<?xml version=': 0, '>': 187, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 54, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6598330', 'n_tokens_mistral': 17140, 'n_tokens_neox': 14508, 'n_words': 3619}