embedding.C File Reference

Go to the source code of this file.

Functions
void	usage_error (const char *progname)
template<typename T >
T	assert_argument (GetPot &cl, const std::string &argname, const char *progname, const T &defaultarg)
int	main (int argc, char **argv)

Function Documentation

assert_argument()

template<typename T >

T assert_argument	(	GetPot &	cl,
		const std::string &	argname,
		const char *	progname,
		const T &	defaultarg
	)

Definition at line 55 of file embedding.C.

References libMesh::err, and usage_error().

Referenced by main().

{
  if (!cl.search(argname))
    {
      libMesh::err << ("No " + argname + " argument found!") << std::endl;
      usage_error(progname);
    }
  return cl.next(defaultarg);
}

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 70 of file embedding.C.

References assert_argument(), libMesh::Elem::build(), libMesh::Elem::embedding_matrix(), libMesh::ReferenceElem::get(), libMesh::TriangleWrapper::init(), int, libMesh::make_range(), libMesh::FEInterface::n_dofs(), n_nodes, libMesh::Elem::parent_bracketing_nodes(), libMesh::Elem::point(), libMesh::Real, libMesh::FEInterface::shape(), and libMesh::TOLERANCE.

{
  LibMeshInit init(argc, argv);

  GetPot cl(argc, argv);

  const std::string elem_type_string =
    assert_argument(cl, "--elem", argv[0], std::string(""));

  const int childnum =
    assert_argument(cl, "--childnum", argv[0], 0);

  const int denominator =
    assert_argument(cl, "--denominator", argv[0], 1);

  const bool diff = cl.search("--diff");

  const ElemType elem_type =
    Utility::string_to_enum<ElemType>(elem_type_string);

  // Getting an embedding matrix isn't a static function, thanks to
  // situations like Tet diagonal selection
  std::unique_ptr<Elem> elem = Elem::build(elem_type);

  const Elem & ref = ReferenceElem::get(elem_type);

  // Lagrange FE for nodal calculations
  FEType fe_type(elem->default_order());

  const unsigned int n_nodes = FEInterface::n_dofs(fe_type, elem.get());
  libmesh_error_msg_if(n_nodes != elem->n_nodes(), "Bad FEInterface value?");

  std::vector<Node> nodes(n_nodes);

  // Get the child vertex positions from childnum; those are easy.
  for (auto v : make_range(elem->n_vertices()))
    {
      for (auto n : make_range(n_nodes))
        {
          const Real embed = ref.embedding_matrix(childnum, v, n);
          if (embed == 1.0)
            {
              nodes[v] = ref.point(n);
              elem->set_node(v, &nodes[v]);
            }
          else if (embed != 0.0)
            libmesh_error_msg("Found fractional embedding on vertex!?");
        }
    }

  // Now figure out the others
  for (auto n : make_range(elem->n_vertices(), n_nodes))
    {
      const auto & pbns = ref.parent_bracketing_nodes(childnum, n);
      if (pbns.empty())
        libmesh_error();

      for (auto pbn : pbns)
        nodes[n] += (ref.point(pbn.first) + ref.point(pbn.second))/2;
      nodes[n] /= pbns.size();
      elem->set_node(n, &nodes[n]);
    }

  const unsigned int denomdigits = std::ceil(std::log10(denominator));
  const unsigned int spacing = denomdigits*2+3;

  std::cout.precision(17);

  std::cout << "    // embedding matrix for child " << childnum << '\n';
  std::cout << "    {\n";
  std::cout << "      //";
  for (auto i : make_range(n_nodes))
    {
      const unsigned int indexdigits =
        std::ceil(std::log10(i));
      const int padding = spacing-indexdigits-2*(i==0);
      for (int k=0; k < padding; ++k)
        std::cout << ' ';
      std::cout << i;
      if (i+1 == n_nodes)
        std::cout << '\n';
      else
        std::cout << ',';
    }

  for (auto i : make_range(n_nodes))
    {
      const Point & pt = elem->point(i);
      std::cout << "      {";
      for (auto j : make_range(n_nodes))
        {
          Real shape = FEInterface::shape(fe_type, elem.get(), j, pt);

          // Don't print -0 or 1e-17; we don't tolerate FP error at 0
          if (std::abs(shape) < TOLERANCE*std::sqrt(TOLERANCE))
            shape = 0;

          const Real embed = ref.embedding_matrix(childnum, i, j);
          if (diff &&
              std::abs(shape - embed) < TOLERANCE*std::sqrt(TOLERANCE))
            {
              for (unsigned int k=0; k != spacing; ++k)
                std::cout << '+';
              if (j+1 != n_nodes)
                std::cout << ',';
            }
          else
            {
              int oldnumerator = int(std::round(shape*denominator));
              int newnumerator = oldnumerator;
              int newdenominator = denominator;

              if (std::abs(oldnumerator-shape*denominator) < TOLERANCE*sqrt(TOLERANCE))
                {
                  int the_gcd = std::gcd(newnumerator, newdenominator);
                  newnumerator /= the_gcd;
                  newdenominator /= the_gcd;
                }

              const unsigned int newdenomdigits =
                std::ceil(std::log10(newdenominator));
              std::ostringstream ostr;
              ostr << (shape*newdenominator);
              const int padding =
                (shape != 0.0 && newdenominator != 1) ?
                int(spacing)-newdenomdigits-2-ostr.str().size() :
                int(spacing)-ostr.str().size();
              for (int k=0; k < padding; ++k)
                std::cout << ' ';
              std::cout << ostr.str();
              if (shape != 0.0 && newdenominator != 1)
                {
                  if (1 << (int)std::round(std::log2(newdenominator)) ==
                      newdenominator)
                    std::cout << "/" << newdenominator << '.';
                  // If we don't have a power of 2 we need to make
                  // sure we're dividing at might-exceed-double Real
                  // precision
                  else
                    std::cout << "/r" << newdenominator;
                }
              if (j+1 != n_nodes)
                std::cout << ',';
            }
        }

      if (i+1 == n_nodes)
        std::cout << "}  ";
      else
        std::cout << "}, ";

      std::cout << "// " << i << '\n';
    }
  std::cout << "    },\n";

  return 0;
}

usage_error()

void usage_error

(

const char *

progname

)

Definition at line 41 of file embedding.C.

References libMesh::out.

Referenced by assert_argument().

{
  libMesh::out << "Options: " << progname << '\n'
               << " --elem type           elem type (e.g. TET14)\n"
               << " --childnum num        child number\n"
               << " --denominator num     denominator to use\n"
               << " --diff                only output diffs from existing\n"
               << std::endl;

  exit(1);
}