diff --git a/examples/vector_fe/vector_fe_ex6/run.sh b/examples/vector_fe/vector_fe_ex6/run.sh index 9d5499a084e..467e29e27fc 100755 --- a/examples/vector_fe/vector_fe_ex6/run.sh +++ b/examples/vector_fe/vector_fe_ex6/run.sh @@ -9,65 +9,65 @@ example_name=vector_fe_ex6 ### Neumann boundary conditions options="dim=2 element_type=TRI6 boundary_condition=neumann" -run_example_no_extra_options "$example_name" "$options" +run_example "$example_name" "$options" options="dim=2 element_type=TRI7 boundary_condition=neumann" -run_example_no_extra_options "$example_name" "$options" +run_example "$example_name" "$options" options="dim=2 element_type=QUAD8 boundary_condition=neumann" -run_example_no_extra_options "$example_name" "$options" +run_example "$example_name" "$options" options="dim=2 element_type=QUAD9 boundary_condition=neumann" -run_example_no_extra_options "$example_name" "$options" +run_example "$example_name" "$options" options="dim=2 order=2 element_type=TRI6 boundary_condition=neumann" -run_example_no_extra_options "$example_name" "$options" +run_example "$example_name" "$options" options="dim=2 order=2 element_type=TRI7 boundary_condition=neumann" -run_example_no_extra_options "$example_name" "$options" +run_example "$example_name" "$options" options="dim=2 order=2 element_type=QUAD8 boundary_condition=neumann" -run_example_no_extra_options "$example_name" "$options" +run_example "$example_name" "$options" options="dim=2 order=2 element_type=QUAD9 boundary_condition=neumann" -run_example_no_extra_options "$example_name" "$options" +run_example "$example_name" "$options" # Subdividing each hex into 24 tets gets expensive in dbg... options="dim=3 element_type=TET14 boundary_condition=neumann grid_size=6" -run_example_no_extra_options "$example_name" "$options" +run_example "$example_name" "$options" options="dim=3 element_type=HEX27 boundary_condition=neumann" -run_example_no_extra_options "$example_name" "$options" +run_example "$example_name" "$options" ### Dirichlet boundary conditions options="dim=2 element_type=TRI6 boundary_condition=dirichlet" -run_example_no_extra_options "$example_name" "$options" +run_example "$example_name" "$options" options="dim=2 element_type=TRI7 boundary_condition=dirichlet" -run_example_no_extra_options "$example_name" "$options" +run_example "$example_name" "$options" options="dim=2 element_type=QUAD8 boundary_condition=dirichlet" -run_example_no_extra_options "$example_name" "$options" +run_example "$example_name" "$options" options="dim=2 element_type=QUAD9 boundary_condition=dirichlet" -run_example_no_extra_options "$example_name" "$options" +run_example "$example_name" "$options" options="dim=2 order=2 element_type=TRI6 boundary_condition=dirichlet" -run_example_no_extra_options "$example_name" "$options" +run_example "$example_name" "$options" options="dim=2 order=2 element_type=TRI7 boundary_condition=dirichlet" -run_example_no_extra_options "$example_name" "$options" +run_example "$example_name" "$options" options="dim=2 order=2 element_type=QUAD8 boundary_condition=dirichlet" -run_example_no_extra_options "$example_name" "$options" +run_example "$example_name" "$options" options="dim=2 order=2 element_type=QUAD9 boundary_condition=dirichlet" -run_example_no_extra_options "$example_name" "$options" +run_example "$example_name" "$options" # Subdividing each hex into 24 tets gets expensive in dbg... options="dim=3 element_type=TET14 boundary_condition=dirichlet grid_size=6" -run_example_no_extra_options "$example_name" "$options" +run_example "$example_name" "$options" options="dim=3 element_type=HEX27 boundary_condition=dirichlet" -run_example_no_extra_options "$example_name" "$options" +run_example "$example_name" "$options" diff --git a/examples/vector_fe/vector_fe_ex9/hdg_problem.C b/examples/vector_fe/vector_fe_ex9/hdg_problem.C index 36461bc6add..bd3f8d7c40f 100644 --- a/examples/vector_fe/vector_fe_ex9/hdg_problem.C +++ b/examples/vector_fe/vector_fe_ex9/hdg_problem.C @@ -175,17 +175,17 @@ HDGProblem::vector_volume_jacobian(DenseMatrix & Jqq, DenseMatrix & u_sol_local, const std::vector & v_sol_local, const unsigned int qp, const unsigned int vel_component) const { - const RealVectorValue U(u_sol_local[qp], v_sol_local[qp]); + const NumberVectorValue U(u_sol_local[qp], v_sol_local[qp]); return U * U(vel_component); } -RealVectorValue +NumberVectorValue HDGProblem::vel_cross_vel_jacobian(const std::vector & u_sol_local, const std::vector & v_sol_local, const unsigned int qp, @@ -194,8 +194,8 @@ HDGProblem::vel_cross_vel_jacobian(const std::vector & u_sol_local, const std::vector> & phi, const unsigned int j) const { - const RealVectorValue U(u_sol_local[qp], v_sol_local[qp]); - RealVectorValue vector_phi_local; + const NumberVectorValue U(u_sol_local[qp], v_sol_local[qp]); + NumberVectorValue vector_phi_local; vector_phi_local(vel_j_component) = phi[j][qp]; auto ret = vector_phi_local * U(vel_component); if (vel_component == vel_j_component) @@ -318,11 +318,11 @@ HDGProblem::pressure_volume_jacobian(DenseMatrix & Jpu, for (const auto j : make_range(scalar_n_dofs)) { { - const Gradient phi((*scalar_phi)[j][qp], 0); + const Gradient phi((*scalar_phi)[j][qp], Number(0)); Jpu(i, j) -= (*JxW)[qp] * ((*grad_scalar_phi)[i][qp] * phi); } { - const Gradient phi(0, (*scalar_phi)[j][qp]); + const Gradient phi(Number(0), (*scalar_phi)[j][qp]); Jpv(i, j) -= (*JxW)[qp] * ((*grad_scalar_phi)[i][qp] * phi); } } @@ -359,11 +359,11 @@ HDGProblem::pressure_face_jacobian(DenseMatrix & Jplm_u, DenseMatrix & Jqq, DenseMatrix & Jqs); - RealVectorValue vel_cross_vel_residual(const std::vector & u_sol_local, - const std::vector & v_sol_local, - const unsigned int qp, - const unsigned int vel_component) const; - - RealVectorValue vel_cross_vel_jacobian(const std::vector & u_sol_local, - const std::vector & v_sol_local, - const unsigned int qp, - const unsigned int vel_component, - const unsigned int vel_j_component, - const std::vector> & phi, - const unsigned int j) const; + NumberVectorValue vel_cross_vel_residual(const std::vector & u_sol_local, + const std::vector & v_sol_local, + const unsigned int qp, + const unsigned int vel_component) const; + + NumberVectorValue vel_cross_vel_jacobian(const std::vector & u_sol_local, + const std::vector & v_sol_local, + const unsigned int qp, + const unsigned int vel_component, + const unsigned int vel_j_component, + const std::vector> & phi, + const unsigned int j) const; void scalar_volume_residual(const std::vector & vel_gradient, const unsigned int vel_component, diff --git a/include/Makefile.in b/include/Makefile.in index 1a60f9a65c5..25df2f3dd70 100644 --- a/include/Makefile.in +++ b/include/Makefile.in @@ -857,6 +857,7 @@ include_HEADERS = \ numerics/eigen_preconditioner.h \ numerics/eigen_sparse_matrix.h \ numerics/eigen_sparse_vector.h \ + numerics/fdm_gradient.h \ numerics/fem_function_base.h \ numerics/function_base.h \ numerics/lumped_mass_matrix.h \ diff --git a/include/error_estimation/exact_solution.h b/include/error_estimation/exact_solution.h index 36e6fde45ad..edbbc66dffc 100644 --- a/include/error_estimation/exact_solution.h +++ b/include/error_estimation/exact_solution.h @@ -38,6 +38,7 @@ class EquationSystems; class Parameters; class Mesh; template class FunctionBase; +template class FEMFunctionBase; enum FEMNormType : int; // Is there any way to simplify this? @@ -112,6 +113,12 @@ class ExactSolution */ void attach_exact_values (const std::vector *> & f); + /** + * Clone and attach arbitrary functors which compute the exact + * values of the EquationSystems' solutions at any point. + */ + void attach_exact_values (const std::vector *> & f); + /** * Clone and attach an arbitrary functor which computes the exact * value of the system \p sys_num solution at any point. @@ -119,6 +126,13 @@ class ExactSolution void attach_exact_value (unsigned int sys_num, FunctionBase * f); + /** + * Clone and attach an arbitrary functor which computes the exact + * value of the system \p sys_num solution at any point. + */ + void attach_exact_value (unsigned int sys_num, + FEMFunctionBase * f); + /** * Attach an arbitrary function which computes the exact value of * the solution at any point. @@ -135,6 +149,12 @@ class ExactSolution */ void attach_exact_derivs (const std::vector *> & g); + /** + * Clone and attach arbitrary functors which compute the exact + * gradients of the EquationSystems' solutions at any point. + */ + void attach_exact_derivs (const std::vector *> & g); + /** * Clone and attach an arbitrary functor which computes the exact * gradient of the system \p sys_num solution at any point. @@ -142,6 +162,13 @@ class ExactSolution void attach_exact_deriv (unsigned int sys_num, FunctionBase * g); + /** + * Clone and attach an arbitrary functor which computes the exact + * gradient of the system \p sys_num solution at any point. + */ + void attach_exact_deriv (unsigned int sys_num, + FEMFunctionBase * g); + /** * Attach an arbitrary function which computes the exact gradient of * the solution at any point. @@ -157,6 +184,11 @@ class ExactSolution * second derivatives of the EquationSystems' solutions at any point. */ void attach_exact_hessians (std::vector *> h); + /** + * Clone and attach arbitrary functors which compute the exact + * second derivatives of the EquationSystems' solutions at any point. + */ + void attach_exact_hessians (std::vector *> h); /** * Clone and attach an arbitrary functor which computes the exact @@ -165,6 +197,13 @@ class ExactSolution void attach_exact_hessian (unsigned int sys_num, FunctionBase * h); + /** + * Clone and attach an arbitrary functor which computes the exact + * second derivatives of the system \p sys_num solution at any point. + */ + void attach_exact_hessian (unsigned int sys_num, + FEMFunctionBase * h); + /** * Attach an arbitrary function which computes the exact second * derivatives of the solution at any point. @@ -315,19 +354,19 @@ class ExactSolution * User-provided functors which compute the exact value of the * solution for each system. */ - std::vector>> _exact_values; + std::vector>> _exact_values; /** * User-provided functors which compute the exact derivative of the * solution for each system. */ - std::vector>> _exact_derivs; + std::vector>> _exact_derivs; /** * User-provided functors which compute the exact hessians of the * solution for each system. */ - std::vector>> _exact_hessians; + std::vector>> _exact_hessians; /** * Data structure which stores the errors: diff --git a/include/include_HEADERS b/include/include_HEADERS index c62e53d1062..0c3888479f4 100644 --- a/include/include_HEADERS +++ b/include/include_HEADERS @@ -257,6 +257,7 @@ include_HEADERS = \ numerics/eigen_preconditioner.h \ numerics/eigen_sparse_matrix.h \ numerics/eigen_sparse_vector.h \ + numerics/fdm_gradient.h \ numerics/fem_function_base.h \ numerics/function_base.h \ numerics/lumped_mass_matrix.h \ diff --git a/include/libmesh/Makefile.am b/include/libmesh/Makefile.am index 4a5ad324c0b..d553ee6eef3 100644 --- a/include/libmesh/Makefile.am +++ b/include/libmesh/Makefile.am @@ -248,6 +248,7 @@ BUILT_SOURCES = \ eigen_preconditioner.h \ eigen_sparse_matrix.h \ eigen_sparse_vector.h \ + fdm_gradient.h \ fem_function_base.h \ function_base.h \ laspack_matrix.h \ @@ -1333,6 +1334,9 @@ eigen_sparse_matrix.h: $(top_srcdir)/include/numerics/eigen_sparse_matrix.h eigen_sparse_vector.h: $(top_srcdir)/include/numerics/eigen_sparse_vector.h $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ +fdm_gradient.h: $(top_srcdir)/include/numerics/fdm_gradient.h + $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ + fem_function_base.h: $(top_srcdir)/include/numerics/fem_function_base.h $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ diff --git a/include/libmesh/Makefile.in b/include/libmesh/Makefile.in index 5a5e26aec00..430109db16b 100644 --- a/include/libmesh/Makefile.in +++ b/include/libmesh/Makefile.in @@ -597,9 +597,9 @@ BUILT_SOURCES = auto_ptr.h dirichlet_boundaries.h dof_map.h \ dense_subvector.h dense_vector.h dense_vector_base.h \ diagonal_matrix.h distributed_vector.h eigen_core_support.h \ eigen_preconditioner.h eigen_sparse_matrix.h \ - eigen_sparse_vector.h fem_function_base.h function_base.h \ - laspack_matrix.h laspack_vector.h lumped_mass_matrix.h \ - numeric_vector.h parsed_fem_function.h \ + eigen_sparse_vector.h fdm_gradient.h fem_function_base.h \ + function_base.h laspack_matrix.h laspack_vector.h \ + lumped_mass_matrix.h numeric_vector.h parsed_fem_function.h \ parsed_fem_function_parameter.h parsed_function.h \ parsed_function_parameter.h petsc_macro.h petsc_matrix.h \ petsc_matrix_base.h petsc_matrix_shell_matrix.h \ @@ -1666,6 +1666,9 @@ eigen_sparse_matrix.h: $(top_srcdir)/include/numerics/eigen_sparse_matrix.h eigen_sparse_vector.h: $(top_srcdir)/include/numerics/eigen_sparse_vector.h $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ +fdm_gradient.h: $(top_srcdir)/include/numerics/fdm_gradient.h + $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ + fem_function_base.h: $(top_srcdir)/include/numerics/fem_function_base.h $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ diff --git a/include/numerics/fdm_gradient.h b/include/numerics/fdm_gradient.h new file mode 100644 index 00000000000..23439ba74c5 --- /dev/null +++ b/include/numerics/fdm_gradient.h @@ -0,0 +1,153 @@ +// The libMesh Finite Element Library. +// Copyright (C) 2002-2024 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner + +// 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 + +#ifndef LIBMESH_FDM_GRADIENT_H +#define LIBMESH_FDM_GRADIENT_H + +// libMesh includes +#include "libmesh/fem_function_base.h" +#include "libmesh/libmesh_common.h" +#include "libmesh/tensor_tools.h" + +namespace libMesh +{ + +template +class FDMGradient : public FEMFunctionBase +{ +public: + typedef typename TensorTools::DecrementRank::type ValType; + + FDMGradient (FEMFunctionBase & value_func, + Real eps) : + _val_func(value_func.clone()), _eps(eps) + {} + + virtual void init_context (const FEMContext & c) override + { _val_func->init_context(c); } + + virtual std::unique_ptr> clone () const override + { return std::make_unique>(*_val_func, _eps); } + + virtual GradType operator() (const FEMContext & c, + const Point & p, + const Real time = 0.) override + { + GradType g; + + auto & val = *_val_func; + + Real one_over_dim = Real(0.5) / _eps; + + g(0) = (val(c, p+Point(_eps), time) - + val(c, p+Point(-_eps), time)) * one_over_dim; +#if LIBMESH_DIM > 1 + g(1) = (val(c, p+Point(0,_eps), time) - + val(c, p+Point(0,-_eps), time)) * one_over_dim; +#endif +#if LIBMESH_DIM > 2 + g(2) = (val(c, p+Point(0,0,_eps), time) - + val(c, p+Point(0,0,-_eps), time)) * one_over_dim; +#endif + + return g; + } + + virtual void operator() (const FEMContext & c, + const Point & p, + const Real time, + DenseVector & output) override + { + auto sz = output.size(); + DenseVector v(sz); + + auto & val = *_val_func; + + val(c, p+Point(_eps), time, v); + for (auto i : make_range(sz)) + output(i)(0) = v(i); + + val(c, p+Point(-_eps), time, v); + for (auto i : make_range(sz)) + { + output(i)(0) -= v(i); + output(i)(0) /= 2; + } + +#if LIBMESH_DIM > 1 + val(c, p+Point(0,_eps), time, v); + for (auto i : make_range(sz)) + output(i)(1) = v(i); + + val(c, p+Point(0,-_eps), time, v); + for (auto i : make_range(sz)) + { + output(i)(1) -= v(i); + output(i)(1) /= 2; + } +#endif +#if LIBMESH_DIM > 2 + val(c, p+Point(0,0,_eps), time, v); + for (auto i : make_range(sz)) + output(i)(2) = v(i); + + val(c, p+Point(0,0,-_eps), time, v); + for (auto i : make_range(sz)) + { + output(i)(2) -= v(i); + output(i)(2) /= 2; + } +#endif + } + + + virtual GradType component (const FEMContext & c, + unsigned int i, + const Point & p, + Real time) override + { + GradType g; + + auto & val = *_val_func; + + Real one_over_dim = Real(0.5) / _eps; + + g(0) = (val.component(c, i, p+Point(_eps), time) - + val.component(c, i, p+Point(-_eps), time)) * one_over_dim; +#if LIBMESH_DIM > 1 + g(1) = (val.component(c, i, p+Point(0,_eps), time) - + val.component(c, i, p+Point(0,-_eps), time)) * one_over_dim; +#endif +#if LIBMESH_DIM > 2 + g(2) = (val.component(c, i, p+Point(0,0,_eps), time) - + val.component(c, i, p+Point(0,0,-_eps), time)) * one_over_dim; +#endif + + return g; + } + +private: + + std::unique_ptr> _val_func; + + Real _eps; +}; + + +} // namespace libMesh + +#endif // LIBMESH_FDM_GRADIENT_H diff --git a/include/numerics/fem_function_base.h b/include/numerics/fem_function_base.h index 88031749697..bd577a22019 100644 --- a/include/numerics/fem_function_base.h +++ b/include/numerics/fem_function_base.h @@ -63,6 +63,11 @@ class FEMFunctionBase FEMFunctionBase & operator= (FEMFunctionBase &&) = default; virtual ~FEMFunctionBase () = default; + /** + * Any post-construction initialization + */ + virtual void init () {} + /** * Prepares a context object for use. * diff --git a/include/numerics/parsed_fem_function.h b/include/numerics/parsed_fem_function.h index d1b809634e9..d33001b53df 100644 --- a/include/numerics/parsed_fem_function.h +++ b/include/numerics/parsed_fem_function.h @@ -396,14 +396,26 @@ ParsedFEMFunction::init_context (const FEMContext & c) { FEBase * elem_fe; c.get_element_fe(v, elem_fe); + bool request_nothing = true; if (_n_requested_vars) - elem_fe->get_phi(); + { + elem_fe->get_phi(); + request_nothing = false; + } if (_n_requested_grad_components) - elem_fe->get_dphi(); + { + elem_fe->get_dphi(); + request_nothing = false; + } #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES if (_n_requested_hess_components) - elem_fe->get_d2phi(); + { + elem_fe->get_d2phi(); + request_nothing = false; + } #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES + if (request_nothing) + elem_fe->get_nothing(); } if (_requested_normals) diff --git a/include/numerics/wrapped_functor.h b/include/numerics/wrapped_functor.h index e17748aa6ce..7ac481417fa 100644 --- a/include/numerics/wrapped_functor.h +++ b/include/numerics/wrapped_functor.h @@ -23,7 +23,9 @@ // Local Includes #include "libmesh/fem_function_base.h" #include "libmesh/function_base.h" +#include "libmesh/int_range.h" #include "libmesh/point.h" +#include "libmesh/wrapped_function.h" // C++ includes #include @@ -54,6 +56,18 @@ class WrappedFunctor : public FEMFunctionBase : _func(func.clone()) { } + /** + * Constructor to wrap scalar-valued function pointers. + */ + WrappedFunctor (const System & sys, + Output fptr(const Point & p, + const Parameters & parameters, + const std::string & sys_name, + const std::string & unknown_name) = nullptr, + const Parameters * parameters = nullptr, + unsigned int varnum=0) : + _func(std::make_unique>(sys, fptr, parameters, varnum)) {} + /** * This class can't be copy constructed or assigned because it * contains a unique_ptr member. @@ -68,6 +82,16 @@ class WrappedFunctor : public FEMFunctionBase WrappedFunctor & operator= (WrappedFunctor &&) = default; virtual ~WrappedFunctor () = default; + /** + * Any post-construction initialization + */ + virtual void init () override { _func->init(); } + + /** + * Tell the context we don't need anything from it + */ + virtual void init_context (const FEMContext & c) override; + virtual std::unique_ptr> clone () const override { return std::make_unique>(*_func); @@ -96,6 +120,20 @@ class WrappedFunctor : public FEMFunctionBase }; +template +void WrappedFunctor::init_context (const FEMContext & c) +{ + for (auto dim : c.elem_dimensions()) + { + for (auto v : make_range(c.n_vars())) + { + FEAbstract * fe; + c.get_element_fe(v, fe, dim); + fe->get_nothing(); + } + } +} + } // namespace libMesh diff --git a/src/apps/L2system.C b/src/apps/L2system.C index 4e95388203c..f082d317eac 100644 --- a/src/apps/L2system.C +++ b/src/apps/L2system.C @@ -29,9 +29,9 @@ using namespace libMesh; -L2System::~L2System () = default; +HilbertSystem::~HilbertSystem () = default; -void L2System::init_data () +void HilbertSystem::init_data () { this->add_variable ("u", static_cast(_fe_order), Utility::string_to_enum(_fe_family)); @@ -42,7 +42,7 @@ void L2System::init_data () -void L2System::init_context(DiffContext & context) +void HilbertSystem::init_context(DiffContext & context) { FEMContext & c = cast_ref(context); @@ -63,6 +63,9 @@ void L2System::init_context(DiffContext & context) my_fe->get_phi(); my_fe->get_xyz(); + if (this->_hilbert_order > 0) + my_fe->get_dphi(); + c.get_side_fe( 0, my_fe, dim ); my_fe->get_nothing(); } @@ -74,16 +77,16 @@ void L2System::init_context(DiffContext & context) { input_context = std::make_unique(*input_system); - libmesh_assert(goal_func.get()); - goal_func->init_context(*input_context); + libmesh_assert(_goal_func.get()); + _goal_func->init_context(*input_context); } FEMSystem::init_context(context); } -bool L2System::element_time_derivative (bool request_jacobian, - DiffContext & context) +bool HilbertSystem::element_time_derivative (bool request_jacobian, + DiffContext & context) { FEMContext & c = cast_ref(context); @@ -121,12 +124,26 @@ bool L2System::element_time_derivative (bool request_jacobian, for (unsigned int qp=0; qp != n_qpoints; qp++) { - Number u = c.interior_value(0, qp); - - Number ufunc = (*goal_func)(input_c, xyz[qp]); + const Number u = c.interior_value(0, qp); + const Number ufunc = (*_goal_func)(input_c, xyz[qp]); + const Number err_u = u - ufunc; for (unsigned int i=0; i != n_u_dofs; i++) - F(i) += JxW[qp] * ((u - ufunc) * phi[i][qp]); + F(i) += JxW[qp] * (err_u * phi[i][qp]); + + if (_hilbert_order > 0) + { + const std::vector> & dphi = + c.get_element_fe(0)->get_dphi(); + + const Gradient grad_u = c.interior_gradient(0, qp); + Gradient ufuncgrad = (*_goal_grad)(input_c, xyz[qp]); + const Gradient err_grad_u = grad_u - ufuncgrad; + + for (unsigned int i=0; i != n_u_dofs; i++) + F(i) += JxW[qp] * (err_grad_u * dphi[i][qp]); + } + if (request_jacobian) { const Number JxWxD = JxW[qp] * @@ -135,6 +152,16 @@ bool L2System::element_time_derivative (bool request_jacobian, for (unsigned int i=0; i != n_u_dofs; i++) for (unsigned int j=0; j != n_u_dofs; ++j) K(i,j) += JxWxD * (phi[i][qp] * phi[j][qp]); + + if (_hilbert_order > 0) + { + const std::vector> & dphi = + c.get_element_fe(0)->get_dphi(); + + for (unsigned int i=0; i != n_u_dofs; i++) + for (unsigned int j=0; j != n_u_dofs; ++j) + K(i,j) += JxWxD * (dphi[i][qp] * dphi[j][qp]); + } } } // end of the quadrature point qp-loop diff --git a/src/apps/L2system.h b/src/apps/L2system.h index 5e82f6a7381..91b5095844e 100644 --- a/src/apps/L2system.h +++ b/src/apps/L2system.h @@ -17,48 +17,65 @@ // libMesh includes #include "libmesh/enum_fe_family.h" +#include "libmesh/fdm_gradient.h" #include "libmesh/fem_function_base.h" #include "libmesh/fem_system.h" +#include "libmesh/libmesh_common.h" // C++ includes #include #include + // FEMSystem, TimeSolver and NewtonSolver will handle most tasks, // but we must specify element residuals -class L2System : public libMesh::FEMSystem +class HilbertSystem : public libMesh::FEMSystem { public: // Constructor - L2System(libMesh::EquationSystems & es, - const std::string & name, - const unsigned int number) + HilbertSystem(libMesh::EquationSystems & es, + const std::string & name, + const unsigned int number) : libMesh::FEMSystem(es, name, number), input_system(nullptr), _fe_family("LAGRANGE"), _fe_order(1), + _hilbert_order(0), _subdomains_list() {} - // Destructor; deletes extra context objects - ~L2System(); + // Default destructor + ~HilbertSystem(); - std::string & fe_family() { return _fe_family; } - unsigned int & fe_order() { return _fe_order; } + std::string & fe_family() { return _fe_family; } + unsigned int & fe_order() { return _fe_order; } std::set & subdomains_list() { return _subdomains_list; } + unsigned int & hilbert_order() { return _hilbert_order; } + void set_fdm_eps(libMesh::Real eps) { + _fdm_eps = eps; + if (_goal_func.get()) + _goal_grad = std::make_unique>(*_goal_func, _fdm_eps); + } + + void set_goal_func(libMesh::FEMFunctionBase & goal) + { + _goal_func = goal.clone(); + _goal_grad = std::make_unique>(*_goal_func, _fdm_eps); + } + // We want to be able to project functions based on *other* systems' // values. For that we need not only a FEMFunction but also a // reference to the system where it applies and a separate context // object (or multiple separate context objects, in the threaded // case) for that system. - std::unique_ptr > goal_func; - libMesh::System * input_system; +protected: + std::unique_ptr > _goal_func; + std::map> input_contexts; -protected: // System initialization virtual void init_data (); @@ -74,6 +91,17 @@ class L2System : public libMesh::FEMSystem std::string _fe_family; unsigned int _fe_order; + // The Hilbert order our subclass will project with + unsigned int _hilbert_order; + + // The function we will call to finite difference our goal + // function + std::unique_ptr> _goal_grad; + + // The perturbation we will use when finite differencing our goal + // function + libMesh::Real _fdm_eps; + // Which subdomains to integrate on (all subdomains, if empty()) std::set _subdomains_list; }; diff --git a/src/apps/calculator.C b/src/apps/calculator.C index c8e652d2b44..2e46198fdd8 100644 --- a/src/apps/calculator.C +++ b/src/apps/calculator.C @@ -18,14 +18,16 @@ // Open the input mesh and corresponding solution file named in command line // arguments, parse the function specified in a command line argument, -// L2-project its value onto the mesh, and output the new solution. +// Project its value onto the mesh, and output the new solution. // libMesh includes #include "L2system.h" #include "libmesh/libmesh.h" #include "libmesh/dof_map.h" +#include "libmesh/enum_norm_type.h" #include "libmesh/equation_systems.h" +#include "libmesh/exact_solution.h" #include "libmesh/getpot.h" #include "libmesh/mesh.h" #include "libmesh/newton_solver.h" @@ -33,6 +35,7 @@ #include "libmesh/parsed_fem_function.h" #include "libmesh/parsed_function.h" #include "libmesh/point.h" +#include "libmesh/sparse_matrix.h" #include "libmesh/steady_solver.h" #include "libmesh/wrapped_functor.h" #include "libmesh/elem.h" @@ -52,16 +55,21 @@ using namespace libMesh; void usage_error(const char * progname) { libMesh::out << "Options: " << progname << '\n' - << " --dim d mesh dimension [default: autodetect]\n" + << " --dim d mesh dimension [default: autodetect]\n" << " --inmesh filename input mesh file\n" + << " --inmat filename input constraint matrix [default: none]\n" << " --insoln filename input solution file\n" << " --calc func function to calculate\n" - << " --insys num input system number [default: 0]\n" - << " --outsoln filename output solution file [default: out_]\n" - << " --family famname output FEM family [default: LAGRANGE]\n" - << " --order num output FEM order [default: 1]\n" - << " --subdomain num subdomain id restriction [default: all subdomains]\n" - << " --integral only calculate integral, not projection\n" + << " --insys sysnum input system number [default: 0]\n" + << " --outsoln filename output solution file [default: out_]\n" + << " --family famname output FEM family [default: LAGRANGE]\n" + << " --order p output FEM order [default: 1]\n" + << " --subdomain sbd_id each subdomain to check [default: all subdomains]\n" + << " --hilbert order Hilbert space to use [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" + << " --integral only calculate func integral, not projection\n" << std::endl; exit(1); @@ -161,15 +169,31 @@ int main(int argc, char ** argv) const std::string meshname = assert_argument(cl, "--inmesh", argv[0], std::string("")); - const std::string solnname = cl.follow(std::string(""), "--insoln"); + libMesh::out << "Reading mesh " << meshname << std::endl; + old_mesh.read(meshname); + + const std::string matname = + cl.follow(std::string(""), "--inmat"); - if (meshname != "") + if (matname != "") { - old_mesh.read(meshname); - std::cout << "Old Mesh:" << std::endl; - old_mesh.print_info(); + 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::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(); + + const std::string solnname = cl.follow(std::string(""), "--insoln"); + // Load the old solution from --insoln filename, if that's been // specified. EquationSystems old_es(old_mesh); @@ -181,12 +205,14 @@ int main(int argc, char ** argv) const std::string family = cl.follow(std::string("LAGRANGE"), "--family"); - const int order = cl.follow(1, "--order"); + const unsigned int order = cl.follow(1u, "--order"); std::unique_ptr> goal_function; if (solnname != "") { + libMesh::out << "Reading solution " << solnname << std::endl; + old_es.read(solnname, EquationSystems::READ_HEADER | EquationSystems::READ_DATA | @@ -245,6 +271,9 @@ int main(int argc, char ** argv) const std::string outsolnname = cl.follow(default_outsolnname, "--outsoln"); + // Output results in high precision + libMesh::out << std::setprecision(std::numeric_limits::max_digits10); + if (!cl.search("--integral")) { // Create a new mesh and a new EquationSystems @@ -253,7 +282,9 @@ int main(int argc, char ** argv) EquationSystems new_es(new_mesh); - L2System & new_sys = new_es.add_system(current_sys_name); + HilbertSystem & new_sys = new_es.add_system(current_sys_name); + + new_sys.hilbert_order() = cl.follow(0u, "--hilbert"); new_sys.subdomains_list() = std::move(subdomains_list); @@ -263,7 +294,11 @@ int main(int argc, char ** argv) new_sys.fe_family() = family; new_sys.fe_order() = order; - new_sys.goal_func = goal_function->clone(); + new_sys.set_goal_func(*goal_function); + + Real fdm_eps = cl.follow(Real(TOLERANCE), "--fdm_eps"); + + new_sys.set_fdm_eps(fdm_eps); if (solnname != "") new_sys.input_system = &old_es.get_system(current_sys_name); @@ -277,6 +312,44 @@ int main(int argc, char ** argv) new_sys.solve(); + // Calculate the error if requested + if (cl.search(1, "--error_q")) + { + const unsigned int error_q = cl.next(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 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 << "L2 error in " << var_name << ": " << + exact_sol.h1_error(current_sys_name, var_name) << + std::endl; + } + } + } + // Write out the new solution file new_es.write(outsolnname.c_str(), EquationSystems::WRITE_DATA | @@ -299,7 +372,7 @@ int main(int argc, char ** argv) Number integral = integrate.integral(); old_mesh.comm().sum(integral); - std::cout << "Integral is " << integral << std::endl; + libMesh::out << "Integral is " << integral << std::endl; } return 0; diff --git a/src/error_estimation/exact_solution.C b/src/error_estimation/exact_solution.C index 468f30f6e74..5b8b53df073 100644 --- a/src/error_estimation/exact_solution.C +++ b/src/error_estimation/exact_solution.C @@ -29,6 +29,7 @@ #include "libmesh/parallel.h" #include "libmesh/quadrature.h" #include "libmesh/wrapped_function.h" +#include "libmesh/wrapped_functor.h" #include "libmesh/fe_interface.h" #include "libmesh/raw_accessor.h" #include "libmesh/tensor_tools.h" @@ -105,7 +106,7 @@ void ExactSolution::attach_exact_value (ValueFunctionPointer fptr) for (auto sys : make_range(_equation_systems.n_systems())) { const System & system = _equation_systems.get_system(sys); - _exact_values.emplace_back(std::make_unique>(system, fptr, &_equation_systems.parameters)); + _exact_values.emplace_back(std::make_unique>(system, fptr, &_equation_systems.parameters)); } // If we're using exact values, we're not using a fine grid solution @@ -119,6 +120,17 @@ void ExactSolution::attach_exact_values (const std::vector // entry for each system. _exact_values.clear(); + for (auto ptr : f) + _exact_values.emplace_back(ptr ? std::make_unique>(*ptr) : nullptr); +} + + +void ExactSolution::attach_exact_values (const std::vector *> & f) +{ + // Automatically delete any previous _exact_values entries, then add a new + // entry for each system. + _exact_values.clear(); + for (auto ptr : f) _exact_values.emplace_back(ptr ? ptr->clone() : nullptr); } @@ -130,6 +142,17 @@ void ExactSolution::attach_exact_value (unsigned int sys_num, if (_exact_values.size() <= sys_num) _exact_values.resize(sys_num+1); + if (f) + _exact_values[sys_num] = std::make_unique>(*f); +} + + +void ExactSolution::attach_exact_value (unsigned int sys_num, + FEMFunctionBase * f) +{ + if (_exact_values.size() <= sys_num) + _exact_values.resize(sys_num+1); + if (f) _exact_values[sys_num] = f->clone(); } @@ -146,7 +169,7 @@ void ExactSolution::attach_exact_deriv (GradientFunctionPointer gptr) for (auto sys : make_range(_equation_systems.n_systems())) { const System & system = _equation_systems.get_system(sys); - _exact_derivs.emplace_back(std::make_unique>(system, gptr, &_equation_systems.parameters)); + _exact_derivs.emplace_back(std::make_unique>(system, gptr, &_equation_systems.parameters)); } // If we're using exact values, we're not using a fine grid solution @@ -160,6 +183,17 @@ void ExactSolution::attach_exact_derivs (const std::vector>(*ptr) : nullptr); +} + + +void ExactSolution::attach_exact_derivs (const std::vector *> & g) +{ + // Automatically delete any previous _exact_derivs entries, then add a new + // entry for each system. + _exact_derivs.clear(); + for (auto ptr : g) _exact_derivs.emplace_back(ptr ? ptr->clone() : nullptr); } @@ -171,6 +205,17 @@ void ExactSolution::attach_exact_deriv (unsigned int sys_num, if (_exact_derivs.size() <= sys_num) _exact_derivs.resize(sys_num+1); + if (g) + _exact_derivs[sys_num] = std::make_unique>(*g); +} + + +void ExactSolution::attach_exact_deriv (unsigned int sys_num, + FEMFunctionBase * g) +{ + if (_exact_derivs.size() <= sys_num) + _exact_derivs.resize(sys_num+1); + if (g) _exact_derivs[sys_num] = g->clone(); } @@ -187,7 +232,7 @@ void ExactSolution::attach_exact_hessian (HessianFunctionPointer hptr) for (auto sys : make_range(_equation_systems.n_systems())) { const System & system = _equation_systems.get_system(sys); - _exact_hessians.emplace_back(std::make_unique>(system, hptr, &_equation_systems.parameters)); + _exact_hessians.emplace_back(std::make_unique>(system, hptr, &_equation_systems.parameters)); } // If we're using exact values, we're not using a fine grid solution @@ -201,6 +246,17 @@ void ExactSolution::attach_exact_hessians (std::vector *> h // entry for each system. _exact_hessians.clear(); + for (auto ptr : h) + _exact_hessians.emplace_back(ptr ? std::make_unique>(*ptr) : nullptr); +} + + +void ExactSolution::attach_exact_hessians (std::vector *> h) +{ + // Automatically delete any previous _exact_hessians entries, then add a new + // entry for each system. + _exact_hessians.clear(); + for (auto ptr : h) _exact_hessians.emplace_back(ptr ? ptr->clone() : nullptr); } @@ -212,6 +268,17 @@ void ExactSolution::attach_exact_hessian (unsigned int sys_num, if (_exact_hessians.size() <= sys_num) _exact_hessians.resize(sys_num+1); + if (h) + _exact_hessians[sys_num] = std::make_unique>(*h); +} + + +void ExactSolution::attach_exact_hessian (unsigned int sys_num, + FEMFunctionBase * h) +{ + if (_exact_hessians.size() <= sys_num) + _exact_hessians.resize(sys_num+1); + if (h) _exact_hessians[sys_num] = h->clone(); } @@ -469,6 +536,8 @@ void ExactSolution::_compute_error(std::string_view sys_name, _equation_systems_fine->get_system(sys_name) : _equation_systems.get_system (sys_name); + FEMContext context(computed_system); + const MeshBase & mesh = computed_system.get_mesh(); const Real time = _equation_systems.get_system(sys_name).time; @@ -509,25 +578,50 @@ void ExactSolution::_compute_error(std::string_view sys_name, coarse_values->init(); } + // Grab which element dimensions are present in the mesh + const std::set & elem_dims = mesh.elem_dimensions(); + // Initialize any functors we're going to use for (auto & ev : _exact_values) if (ev) - ev->init(); + { + ev->init(); + ev->init_context(context); + } for (auto & ed : _exact_derivs) if (ed) - ed->init(); + { + ed->init(); + ed->init_context(context); + } for (auto & eh : _exact_hessians) if (eh) - eh->init(); + { + eh->init(); + eh->init_context(context); + } + + // If we have *no* functors we intend to use (because we're using a + // fine system, or because our exact solution is zero and we're just + // computing norms in an outdated way) then let our FE objects know + // we don't actually need anything from them, so they don't think + // we've just invoked them in a deprecated "compute everything" + // fashion. + if (_exact_values.empty() && _exact_derivs.empty() && + _exact_hessians.empty()) + for (auto dim : elem_dims) + for (auto v : make_range(computed_system.n_vars())) + { + FEAbstract * fe; + context.get_element_fe(v, fe, dim); + fe->get_nothing(); + } // Get a reference to the dofmap and mesh for that system const DofMap & computed_dof_map = computed_system.get_dof_map(); - // Grab which element dimensions are present in the mesh - const std::set & elem_dims = mesh.elem_dimensions(); - // Zero the error before summation // 0 - sum of square of function error (L2) // 1 - sum of square of gradient error (H1 semi) @@ -652,6 +746,9 @@ void ExactSolution::_compute_error(std::string_view sys_name, // for the current element fe->reinit (elem); + context.pre_fe_reinit(computed_system, elem); + context.elem_fe_reinit(); + // Get the local to global degree of freedom maps computed_dof_map.dof_indices (elem, dof_indices, var); @@ -709,7 +806,7 @@ void ExactSolution::_compute_error(std::string_view sys_name, for (unsigned int c = 0; c < n_vec_dim; c++) exact_val_accessor(c) = _exact_values[sys_num]-> - component(var_component+c, q_point[qp], time); + component(context, var_component+c, q_point[qp], time); } else if (_equation_systems_fine) { @@ -739,7 +836,7 @@ void ExactSolution::_compute_error(std::string_view sys_name, for (unsigned int d = 0; d < LIBMESH_DIM; d++) exact_grad_accessor(d + c*LIBMESH_DIM) = _exact_derivs[sys_num]-> - component(var_component+c, q_point[qp], time)(d); + component(context, var_component+c, q_point[qp], time)(d); } else if (_equation_systems_fine) { @@ -809,7 +906,7 @@ void ExactSolution::_compute_error(std::string_view sys_name, for (unsigned int e =0; e < dim; e++) exact_hess_accessor(d + e*dim + c*dim*dim) = _exact_hessians[sys_num]-> - component(var_component+c, q_point[qp], time)(d,e); + component(context, var_component+c, q_point[qp], time)(d,e); // FIXME: operator- is not currently implemented for TypeNTensor const typename FEGenericBase::OutputNumberTensor grad2_error = grad2_u_h - exact_hess; diff --git a/src/fe/fe_abstract.C b/src/fe/fe_abstract.C index f8d06cb4e9b..0f1e07052c1 100644 --- a/src/fe/fe_abstract.C +++ b/src/fe/fe_abstract.C @@ -140,6 +140,18 @@ std::unique_ptr FEAbstract::build(const unsigned int dim, case SCALAR: return std::make_unique>(fet); + case NEDELEC_ONE: + return std::make_unique>(fet); + + case RAVIART_THOMAS: + return std::make_unique>(fet); + + case L2_RAVIART_THOMAS: + return std::make_unique>(fet); + + case SUBDIVISION: + return std::make_unique(fet); + default: libmesh_error_msg("ERROR: Bad FEType.family= " << Utility::enum_to_string(fet.family)); } @@ -205,6 +217,18 @@ std::unique_ptr FEAbstract::build(const unsigned int dim, case SCALAR: return std::make_unique>(fet); + case NEDELEC_ONE: + return std::make_unique>(fet); + + case RAVIART_THOMAS: + return std::make_unique>(fet); + + case L2_RAVIART_THOMAS: + return std::make_unique>(fet); + + case SUBDIVISION: + return std::make_unique(fet); + default: libmesh_error_msg("ERROR: Bad FEType.family= " << Utility::enum_to_string(fet.family)); } diff --git a/src/numerics/sparse_matrix.C b/src/numerics/sparse_matrix.C index 18ef4927e8f..ef3f649794e 100644 --- a/src/numerics/sparse_matrix.C +++ b/src/numerics/sparse_matrix.C @@ -581,6 +581,14 @@ void SparseMatrix::read_matlab(const std::string & filename) std::size_t m = 0, n = 0; + // If we don't already have this size, we'll need to reinit, and + // determine which rows+columns each processor is in charge of. + std::vector new_row_starts, new_row_stops, + new_col_starts, new_col_stops; + + numeric_index_type new_row_start, new_row_stop, + new_col_start, new_col_stop; + // We'll read through the file three times: once to get a reliable // value for the matrix size (so we can divvy it up among // processors), then again to get the sparsity to send to each @@ -601,6 +609,8 @@ void SparseMatrix::read_matlab(const std::string & filename) // even with a plain ifstream. What happened to disk caching!? std::vector> entries; + // First read through the file, saving size and entry data + { // We'll read the matrix on processor 0 rather than try to juggle // parallel I/O. if (this->processor_id() == 0) @@ -609,8 +619,6 @@ void SparseMatrix::read_matlab(const std::string & filename) // more robust to formatting. const std::regex start_regex // assignment like "zzz = [" ("\\s*\\w+\\s*=\\s*\\["); - const std::regex entry_regex // row/col/val like "1 1 -2.0e-4" - ("(\\d+)\\s+(\\d+)\\s+([+-]?(\\d+([.]\\d*)?([eE][+-]?\\d+)?|[.]\\d+([eE][+-]?\\d+)?))"); const std::regex end_regex // end of assignment ("^[^%]*\\]"); @@ -640,6 +648,7 @@ void SparseMatrix::read_matlab(const std::string & filename) // way to get the size is if it ended up in a comment const std::regex size_regex // comment like "% size = 8 8" ("%\\s*[Ss][Ii][Zz][Ee]\\s*=\\s*(\\d+)\\s+(\\d+)"); + const std::string whitespace = " \t"; bool have_started = false; bool have_ended = false; @@ -653,34 +662,23 @@ void SparseMatrix::read_matlab(const std::string & filename) { std::smatch sm; - if (std::regex_search(line, sm, size_regex)) - { - const std::string msize = sm[1]; - const std::string nsize = sm[2]; - m = std::stoull(msize); - n = std::stoull(nsize); - } + // First, try to match an entry. This is the most common + // case so we won't rely on slow std::regex for it. + // stringstream is at least an improvement over that. - if (std::regex_search(line, start_regex)) - have_started = true; + // Look for row/col/val like "1 1 -2.0e-4" - if (std::regex_search(line, sm, entry_regex)) - { - libmesh_error_msg_if - (!have_started, "Confused by premature entries in matrix file " << filename); + std::istringstream l(line); - const std::string istr = sm[1]; - const std::string jstr = sm[2]; - const std::string cstr = sm[3]; + std::size_t i, j; + T value; - std::size_t i = std::stoull(istr); - std::size_t j = std::stoull(jstr); + l >> i >> j >> value; - // Try to be compatible with - // higher-than-double-precision T - std::stringstream ss(cstr); - T value; - ss >> value; + if (!l.fail()) + { + libmesh_error_msg_if + (!have_started, "Confused by premature entries in matrix file " << filename); entries.emplace_back(cast_int(i), cast_int(j), @@ -701,7 +699,18 @@ void SparseMatrix::read_matlab(const std::string & filename) largest_j_seen = std::max(j, largest_j_seen); } - if (std::regex_search(line, end_regex)) + else if (std::regex_search(line, sm, size_regex)) + { + const std::string msize = sm[1]; + const std::string nsize = sm[2]; + m = std::stoull(msize); + n = std::stoull(nsize); + } + + else if (std::regex_search(line, start_regex)) + have_started = true; + + else if (std::regex_search(line, end_regex)) { have_ended = true; break; @@ -741,17 +750,6 @@ void SparseMatrix::read_matlab(const std::string & filename) this->comm().broadcast(n); } - // If we don't already have this size, we'll need to reinit later, - // and we'll need to redetermine which rows each processor is in - // charge of. - - numeric_index_type - new_row_start = this->processor_id() * m / this->n_processors(), - new_row_stop = (this->processor_id()+1) * m / this->n_processors(); - numeric_index_type - new_col_start = this->processor_id() * n / this->n_processors(), - new_col_stop = (this->processor_id()+1) * n / this->n_processors(); - if (this->initialized() && m == this->m() && n == this->n()) @@ -762,17 +760,27 @@ void SparseMatrix::read_matlab(const std::string & filename) new_col_start = this->col_start(), new_col_stop = this->col_stop(); } + else + { + // Determine which rows/columns each processor will be in charge of + new_row_start = this->processor_id() * m / this->n_processors(), + new_row_stop = (this->processor_id()+1) * m / this->n_processors(); - std::vector new_row_starts, new_row_stops, - new_col_starts, new_col_stops; + new_col_start = this->processor_id() * n / this->n_processors(), + new_col_stop = (this->processor_id()+1) * n / this->n_processors(); + } this->comm().gather(0, new_row_start, new_row_starts); this->comm().gather(0, new_row_stop, new_row_stops); this->comm().gather(0, new_col_start, new_col_starts); this->comm().gather(0, new_col_stop, new_col_stops); - // Reread to deduce the sparsity pattern, or at least the maximum - // number of on- and off- diagonal non-zeros per row. + } // Done reading entry data and broadcasting matrix size + + // Calculate the matrix sparsity and initialize it second + { + // Deduce the sparsity pattern, or at least the maximum number of + // on- and off- diagonal non-zeros per row. numeric_index_type on_diagonal_nonzeros =0, off_diagonal_nonzeros =0; @@ -824,9 +832,9 @@ void SparseMatrix::read_matlab(const std::string & filename) new_col_stop-new_col_start, on_diagonal_nonzeros, off_diagonal_nonzeros); + } - // One last reread to set values - // + // Set the matrix values last. // Convert from 1-based to 0-based indexing if (this->processor_id() == 0) for (auto [i, j, value] : entries)