Classes
struct	Integrate

Functions
void	usage_error (const char *progname)

template<typename T >
T	assert_argument (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	(	const std::string &	argname,
		const char *	progname,
		const T &	defaultarg
	)

Definition at line 91 of file calculator.C.

{
  if (!libMesh::on_command_line(argname))
    {
      libMesh::err << ("No " + argname + " argument found!") << std::endl;
      usage_error(progname);
    }
  return libMesh::command_line_next(argname, defaultarg);
}

References libMesh::command_line_next(), libMesh::err, libMesh::on_command_line(), and usage_error().

Referenced by main().

◆ main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 163 of file calculator.C.

{
  LibMeshInit init(argc, argv);
 
  // In case the mesh file doesn't let us auto-infer dimension, we let
  // the user specify it on the command line
  const unsigned char requested_dim =
    cast_int<unsigned char>(libMesh::command_line_next("--dim", 3));
 
  // Load the old mesh from --inmesh filename.
  // Keep it serialized; we don't want elements on the new mesh to be
  // looking for data on old mesh elements that live off-processor.
  Mesh old_mesh(init.comm(), requested_dim);
 
  const std::string meshname =
    assert_argument("--inmesh", argv[0], std::string(""));
 
  libMesh::out << "Reading mesh " << meshname << std::endl;
  old_mesh.read(meshname);
 
  const std::string matname =
    libMesh::command_line_next("--inmat", std::string(""));
 
  if (matname != "")
    {
      libMesh::out << "Reading matrix " << matname << std::endl;
 
      // For extraction matrices Coreform has been experimenting with
      // PETSc solvers which take the transpose of what we expect, so
      // we'll un-transpose here.
      auto matrix = SparseMatrix<Number>::build (old_mesh.comm());
      matrix->read(matname);
      matrix->get_transpose(*matrix);
 
      old_mesh.copy_constraint_rows(*matrix);
    }
 
  libMesh::out << "Mesh:" << std::endl;
  old_mesh.print_info();
 
  // If we're not using a distributed mesh, this is cheap info to add
  if (old_mesh.is_serial_on_zero())
    {
      const Real mat_tol =
        libMesh::command_line_next("--mattol", Real(0));
 
      const dof_id_type n_components =
        MeshTools::n_connected_components(old_mesh, mat_tol);
      libMesh::out << "Mesh has " << n_components << " connected components." << std::endl;
    }
 
  const std::string solnname =
    libMesh::command_line_next("--insoln", std::string(""));
 
  // Load the old solution from --insoln filename, if that's been
  // specified.
  EquationSystems old_es(old_mesh);
  std::string current_sys_name = "new_sys";
 
  const std::string calcfunc =
    assert_argument("--calc", argv[0], std::string(""));
 
  const std::string family =
    libMesh::command_line_next("--family", std::string("LAGRANGE"));
 
  const unsigned int order = libMesh::command_line_next("--order", 1u);
 
  std::unique_ptr<FEMFunctionBase<Number>> goal_function;
 
  if (solnname != "")
    {
      libMesh::out << "Reading solution " << solnname << std::endl;
 
      old_es.read(solnname,
                  EquationSystems::READ_HEADER |
                  EquationSystems::READ_DATA |
                  EquationSystems::READ_ADDITIONAL_DATA |
                  EquationSystems::READ_BASIC_ONLY);
 
      libMesh::out << "Input system:\n\n";
 
      old_es.print_info();
 
      const unsigned int sysnum =
        libMesh::command_line_next("--insys", 0);
 
      libmesh_assert_less(sysnum, old_es.n_systems());
 
      System & old_sys = old_es.get_system(sysnum);
 
      current_sys_name = old_sys.name();
 
      goal_function =
        std::make_unique<ParsedFEMFunction<Number>>(old_sys, calcfunc);
    }
  else
    {
      current_sys_name = "bare_sys";
 
      // FEMContext doesn't like seeing systems with no variables
      System & dummy_sys = old_es.add_system<System>(current_sys_name);
 
      // And if we're using an IsoGeometric Analysis mesh, we'd better
      // make sure the user can specify order and family to make it
      // iso.
      dummy_sys.add_variable("dummy", static_cast<Order>(order),
                             Utility::string_to_enum<FEFamily>(family));
 
      old_es.init();
 
      goal_function =
        std::make_unique<WrappedFunctor<Number>>(ParsedFunction<Number>(calcfunc));
    }
 
  libMesh::out << "Calculating with system " << current_sys_name << std::endl;
 
  // Subdomains to integrate on
  std::vector<libMesh::subdomain_id_type> subdomain_vec;
  libMesh::command_line_vector("--subdomain", subdomain_vec);
  std::set<libMesh::subdomain_id_type> subdomains_list(subdomain_vec.begin(),
                                                       subdomain_vec.end());
 
  std::string default_outsolnname = "out_soln.xda";
  if (solnname != "")
    default_outsolnname = "out_"+solnname;
  const std::string outsolnname =
    libMesh::command_line_next("--outsoln", default_outsolnname);
 
  // Output results in high precision
  libMesh::out << std::setprecision(std::numeric_limits<Real>::max_digits10);
 
  if (!libMesh::on_command_line("--integral"))
    {
      // Create a new mesh and a new EquationSystems.
      // Copying the old mesh is faster than rereading it, and gets us
      // any constraint rows from an input extraction operator matrix.
      Mesh new_mesh(old_mesh);
 
      EquationSystems new_es(new_mesh);
 
      HilbertSystem & new_sys = new_es.add_system<HilbertSystem>(current_sys_name);
 
      if (libMesh::on_command_line("--print_solutions"))
        new_sys.print_solutions = true;
 
      if (libMesh::on_command_line("--print_residuals"))
        new_sys.print_residuals = true;
 
      if (libMesh::on_command_line("--print_jacobians"))
        new_sys.print_jacobians = true;
 
      if (libMesh::on_command_line("--print_element_solutions"))
        new_sys.print_element_solutions = true;
 
      if (libMesh::on_command_line("--print_element_residuals"))
        new_sys.print_element_residuals = true;
 
      if (libMesh::on_command_line("--print_element_jacobians"))
        new_sys.print_element_jacobians = true;
 
      new_sys.hilbert_order() = libMesh::command_line_next("--hilbert", 0u);
 
      new_sys.subdomains_list() = std::move(subdomains_list);
 
      new_sys.time_solver =
        std::make_unique<SteadySolver>(new_sys);
 
      new_sys.fe_family() = family;
      new_sys.fe_order() = order;
 
      new_sys.set_goal_func(*goal_function);
 
      const Real fdm_eps = libMesh::command_line_next("--fdm_eps", Real(TOLERANCE));
 
      new_sys.set_fdm_eps(fdm_eps);
 
      if (solnname != "")
        new_sys.input_system = &old_es.get_system(current_sys_name);
 
      new_es.init();
 
      libMesh::out << "New system:\n\n";
 
      new_es.print_info();
 
      DiffSolver & solver = *(new_sys.time_solver->diff_solver().get());
      solver.quiet = false;
      solver.verbose = true;
      solver.relative_step_tolerance = 1e-10;
 
      new_sys.solve();
 
      // Integrate the error if requested
      if (libMesh::on_command_line("--error_q"))
        {
          const unsigned int error_q =
            libMesh::command_line_next("--error_q", 0u);
 
          // We just add "u" now but maybe we'll change that
          for (auto v : make_range(new_sys.n_vars()))
            {
              ExactSolution exact_sol(new_es);
              exact_sol.attach_exact_value(0, goal_function.get());
              FDMGradient<Gradient> fdm_gradient(*goal_function, fdm_eps);
              exact_sol.attach_exact_deriv(0, &fdm_gradient);
              exact_sol.extra_quadrature_order(error_q);
 
              const std::string var_name = new_sys.variable_name(v);
              exact_sol.compute_error(current_sys_name, var_name);
 
              libMesh::out << "L2 norm of " << var_name << ": " <<
                new_sys.calculate_norm(*new_sys.solution, v, L2) <<
                std::endl;
 
              libMesh::out << "L2 error in " << var_name << ": " <<
                  exact_sol.l2_error(current_sys_name, var_name) <<
                  std::endl;
 
              if (new_sys.hilbert_order() > 0)
                {
                  libMesh::out << "H1 norm of " << var_name << ": " <<
                    new_sys.calculate_norm(*new_sys.solution, v, H1) <<
                    std::endl;
 
                  libMesh::out << "H1 error in " << var_name << ": " <<
                      exact_sol.h1_error(current_sys_name, var_name) <<
                      std::endl;
                }
            }
        }
 
      // Calculate the jump error indicator if requested
      if (libMesh::on_command_line("--jump_error"))
        {
          const unsigned int jump_error_hilbert =
            libMesh::command_line_next("--jump_error", 0u);
 
          for (auto v : make_range(new_sys.n_vars()))
            {
              std::unique_ptr<JumpErrorEstimator> error_estimator;
              SystemNorm error_norm;
              error_norm.set_weight(0,0);
              error_norm.set_weight(v+1,0);
              error_norm.set_weight(v,1.0);
 
              if (jump_error_hilbert == 0)
                {
                  error_estimator = std::make_unique<DiscontinuityMeasure>();
                  error_norm.set_type(v, L2);
                }
              else if (jump_error_hilbert == 1)
                {
                  error_estimator = std::make_unique<KellyErrorEstimator>();
                  error_norm.set_type(v, H1_SEMINORM);
                }
              else if (jump_error_hilbert == 2)
                {
                  error_estimator = std::make_unique<LaplacianErrorEstimator>();
                  error_norm.set_type(v, H2_SEMINORM);
                }
              else
                libmesh_not_implemented();
 
              if (libMesh::on_command_line("--jump_slits"))
                error_estimator->integrate_slits = true;
 
              ErrorVector error_per_cell;
              error_estimator->error_norm = error_norm;
              error_estimator->estimate_error(new_sys, error_per_cell);
 
              const std::string var_name = new_sys.variable_name(v);
 
              // I really want to just stop supporting libstdc++-8.
              // const Real total_error = std::reduce(error_per_cell.begin(), error_per_cell.end());
 
              ErrorVectorReal total_error = 0;
              for (auto cell_error : error_per_cell)
                total_error += cell_error;
 
              if (jump_error_hilbert == 0)
                libMesh::out << "L2";
              else
                libMesh::out << "H" << jump_error_hilbert;
 
              libMesh::out << " error estimate for " << var_name << ": " <<
                total_error << std::endl;
            }
        }
 
      // Write out the new solution file
      new_es.write(outsolnname.c_str(),
                   EquationSystems::WRITE_DATA |
                   EquationSystems::WRITE_ADDITIONAL_DATA);
      libMesh::out << "Wrote solution " << outsolnname << std::endl;
    }
  else
    {
      // If we aren't doing a projection, then we just do an integral
      System & old_sys = old_es.get_system<System>(current_sys_name);
 
      ConstElemRange active_local_elem_range
        (old_mesh.active_local_elements_begin(),
         old_mesh.active_local_elements_end());
 
      Integrate integrate(old_sys, *goal_function);
 
      Threads::parallel_reduce (active_local_elem_range,
                                integrate);
 
      Number integral = integrate.integral();
      old_mesh.comm().sum(integral);
      libMesh::out << "Integral is " << integral << std::endl;
    }
 
  return 0;
}

◆ usage_error()

void usage_error ( const char * progname )

Definition at line 61 of file calculator.C.

{
  libMesh::out << "Options: " << progname << '\n'
               << " --dim        d         mesh dimension               [default: autodetect]\n"
               << " --inmesh     filename  input mesh file\n"
               << " --inmat      filename  input constraint matrix      [default: none]\n"
               << " --mattol     filename  constraint tolerance when testing mesh connectivity\n"
               << "                                                     [default: 0]\n"
               << " --insoln     filename  input solution file\n"
               << " --calc       func      function to calculate\n"
               << " --insys      sysnum    input system number          [default: 0]\n"
               << " --outsoln    filename  output solution file         [default: out_<insoln>]\n"
               << " --family     famname   output FEM family            [default: LAGRANGE]\n"
               << " --order      p         output FEM order             [default: 1]\n"
               << " --subdomain  \"sbd_ids\" each subdomain to check      [default: all subdomains]\n"
               << " --hilbert    order     Hilbert space to project in  [default: 0 => H0]\n"
               << " --fdm_eps    eps       Central diff for dfunc/dx    [default: " << TOLERANCE << "]\n"
               << " --error_q    extra_q   integrate projection error, with adjusted\n"
               << "                       (extra) quadrature order      [default: off, suggested: 0]\n"
               << " --jump_error order     calculate jump error indicator, for specified\n"
               << "                       Hilbert order                 [default: off]\n"
               << " --jump_slits           calculate jumps across slits [default: off]\n"
               << " --integral             only calculate func integral, not projection\n"
               << std::endl;
 
  exit(1);
}

References libMesh::out, and libMesh::TOLERANCE.

Referenced by assert_argument().

Classes

Functions

Function Documentation

◆ assert_argument()

◆ main()

◆ usage_error()