1 #include <libmesh/equation_systems.h> 2 #include <libmesh/replicated_mesh.h> 3 #include <libmesh/numeric_vector.h> 4 #include <libmesh/dof_map.h> 5 #include <libmesh/elem.h> 6 #include <libmesh/face_quad4.h> 7 #include <libmesh/ghosting_functor.h> 8 #include <libmesh/coupling_matrix.h> 9 #include <libmesh/quadrature_gauss.h> 10 #include <libmesh/linear_implicit_system.h> 11 #include <libmesh/mesh_refinement.h> 12 #include <libmesh/mesh_modification.h> 13 #include <libmesh/partitioner.h> 14 #include <libmesh/sparse_matrix.h> 44 _mesh(system.get_mesh()),
45 _point_locator(system.get_mesh().sub_point_locator()),
46 _coupling_matrix(nullptr)
50 CPPUNIT_ASSERT_EQUAL(2, static_cast<int>(_mesh.n_subdomains()));
53 std::set<subdomain_id_type> ids;
54 _mesh.subdomain_ids(ids);
55 CPPUNIT_ASSERT( ids.find(1) != ids.end() );
56 CPPUNIT_ASSERT( ids.find(2) != ids.end() );
58 _point_locator->enable_out_of_mesh_mode();
64 { _coupling_matrix = std::move(coupling_matrix); }
66 virtual void operator()
70 map_type & coupled_elements)
override 72 std::unique_ptr<PointLocatorBase> sub_point_locator = _mesh.sub_point_locator();
74 for(
const auto & elem :
as_range(range_begin,range_end) )
76 std::set<subdomain_id_type> allowed_subdomains;
78 if( elem->subdomain_id() == 1 )
80 var = _system.variable_number(
"V");
81 allowed_subdomains.insert(2);
83 else if( elem->subdomain_id() == 2 )
85 var = _system.variable_number(
"U");
86 allowed_subdomains.insert(1);
89 CPPUNIT_FAIL(
"Something bad happpend");
92 FEType fe_type = _system.variable_type(var);
95 fe->attach_quadrature_rule (&qrule);
97 const std::vector<libMesh::Point> & qpoints = fe->get_xyz();
101 for (
const auto & qp : qpoints )
103 const Elem * overlapping_elem = (*sub_point_locator)( qp, &allowed_subdomains );
105 if( overlapping_elem && overlapping_elem->
processor_id() != p )
106 coupled_elements.emplace(overlapping_elem, _coupling_matrix.get());
142 virtual std::unique_ptr<Partitioner>
clone ()
const override 144 return std::make_unique<OverlappingTestPartitioner>(*this);
153 const unsigned int n)
override 157 this->single_partition_range (
mesh.active_elements_begin(),
mesh.active_elements_end());
160 libmesh_assert_greater (n, 0);
167 for (
auto & elem :
mesh.active_element_ptr_range() )
169 elem->processor_id() = e;
178 unsigned int n_sub_two_elems =
std::distance(
mesh.active_subdomain_elements_begin(2),
179 mesh.active_subdomain_elements_end(2) );
181 unsigned int n_parts_sub_two = 0;
182 if( n_sub_two_elems < n/2 )
183 n_parts_sub_two = n_sub_two_elems;
185 n_parts_sub_two = n/2;
187 const unsigned int n_parts_sub_one = n - n_parts_sub_two;
189 const dof_id_type sub_two_blk_size = cast_int<dof_id_type>
191 mesh.active_subdomain_elements_end(2) )/n_parts_sub_two );
193 this->assign_proc_id_subdomain(
mesh, 2, sub_two_blk_size, n_parts_sub_two, 0 );
196 const dof_id_type sub_one_blk_size = cast_int<dof_id_type>
198 mesh.active_subdomain_elements_end(1) )/n_parts_sub_one );
200 this->assign_proc_id_subdomain(
mesh, 1, sub_one_blk_size, n_parts_sub_one, n_parts_sub_two );
208 const unsigned int n_parts,
214 if ((e/blksize) < n_parts)
215 elem->processor_id() = offset + cast_int<processor_id_type>(e/blksize);
217 elem->processor_id() = offset;
236 std::unique_ptr<EquationSystems>
_es;
243 _mesh = std::make_unique<ReplicatedMesh>(*TestCommWorld);
245 _mesh->set_mesh_dimension(2);
247 _mesh->add_point(
Point(0.0,0.0),0 );
248 _mesh->add_point(
Point(1.0,0.0),1 );
249 _mesh->add_point(
Point(1.0,1.0),2 );
250 _mesh->add_point(
Point(0.0,1.0),3 );
256 for (
unsigned int n=0; n<4; n++)
257 elem->
set_node(n, _mesh->node_ptr(n));
260 _mesh->add_point(
Point(1.0,2.0),4 );
261 _mesh->add_point(
Point(0.0,2.0),5 );
267 elem->
set_node(0, _mesh->node_ptr(3));
268 elem->
set_node(1, _mesh->node_ptr(2));
269 elem->
set_node(2, _mesh->node_ptr(4));
270 elem->
set_node(3, _mesh->node_ptr(5));
273 _mesh->add_point(
Point(0.0,0.0),6 );
274 _mesh->add_point(
Point(1.0,0.0),7 );
275 _mesh->add_point(
Point(1.0,2.0),8 );
276 _mesh->add_point(
Point(0.0,2.0),9 );
282 elem->
set_node(0, _mesh->node_ptr(6));
283 elem->
set_node(1, _mesh->node_ptr(7));
284 elem->
set_node(2, _mesh->node_ptr(8));
285 elem->
set_node(3, _mesh->node_ptr(9));
288 _mesh->partitioner() = std::make_unique<OverlappingTestPartitioner>();
290 _mesh->prepare_for_use();
292 #ifdef LIBMESH_ENABLE_AMR 293 if (n_refinements > 0)
298 #endif // LIBMESH_ENABLE_AMR 303 _es = std::make_unique<EquationSystems>(*_mesh);
306 std::set<subdomain_id_type> sub_one;
309 std::set<subdomain_id_type> sub_two;
329 System & system = _es->get_system(
"SimpleSystem");
331 coupling = std::make_unique<CouplingMatrix>(system.
n_vars());
340 (*coupling)(u_var,v_var) =
true;
341 (*coupling)(l_var,v_var) =
true;
342 (*coupling)(l_var,p_var) =
true;
343 (*coupling)(v_var,u_var) =
true;
344 (*coupling)(v_var,l_var) =
true;
359 #ifdef LIBMESH_HAVE_SOLVER 360 CPPUNIT_TEST( checkCouplingFunctorQuad );
361 CPPUNIT_TEST( checkCouplingFunctorQuadUnifRef );
362 CPPUNIT_TEST( checkCouplingFunctorTri );
363 CPPUNIT_TEST( checkCouplingFunctorTriUnifRef );
364 CPPUNIT_TEST( checkOverlappingPartitioner );
365 CPPUNIT_TEST( checkOverlappingPartitionerUnifRef );
368 CPPUNIT_TEST_SUITE_END();
375 this->build_quad_mesh();
383 { this->run_coupling_functor_test(0); }
386 { this->run_coupling_functor_test(1); }
392 this->run_coupling_functor_test(0);
399 this->run_coupling_functor_test(1);
404 this->run_partitioner_test(0);
409 this->run_partitioner_test(1);
420 #if defined(LIBMESH_ENABLE_AMR) && defined(LIBMESH_HAVE_SOLVER) 421 if( n_refinements > 0 )
429 System & system = _es->get_system(
"SimpleSystem");
436 coupling_functor( _mesh->active_subdomain_elements_begin(1),
437 _mesh->active_subdomain_elements_end(1),
439 subdomain_one_couplings );
441 coupling_functor( _mesh->active_subdomain_elements_begin(2),
442 _mesh->active_subdomain_elements_end(2),
444 subdomain_two_couplings );
447 _mesh->active_subdomain_elements_end(2) );
449 CPPUNIT_ASSERT_EQUAL(n_elems_subdomain_two, static_cast<dof_id_type>(subdomain_one_couplings.size()));
450 CPPUNIT_ASSERT_EQUAL(static_cast<dof_id_type>(2*n_elems_subdomain_two), static_cast<dof_id_type>(subdomain_two_couplings.size()));
455 #ifdef LIBMESH_ENABLE_AMR 456 if( n_refinements > 0 )
462 #endif // LIBMESH_ENABLE_AMR 464 System & system = _es->get_system(
"SimpleSystem");
466 std::set<processor_id_type> sub_one_proc_ids, sub_two_proc_ids;
467 for (
auto & elem : _mesh->active_subdomain_elements_ptr_range(1))
468 sub_one_proc_ids.insert(elem->processor_id());
470 for (
auto & elem : _mesh->active_subdomain_elements_ptr_range(2))
471 sub_two_proc_ids.insert(elem->processor_id());
477 CPPUNIT_ASSERT_EQUAL(1, static_cast<int>(sub_one_proc_ids.size()));
478 CPPUNIT_ASSERT_EQUAL(1, static_cast<int>(sub_two_proc_ids.size()));
479 CPPUNIT_ASSERT( sub_one_proc_ids == sub_two_proc_ids );
485 for (
auto &
id : sub_one_proc_ids )
486 CPPUNIT_ASSERT( sub_two_proc_ids.find(
id) == sub_two_proc_ids.end() );
489 for (
auto &
id : sub_two_proc_ids )
490 CPPUNIT_ASSERT( sub_one_proc_ids.find(
id) == sub_one_proc_ids.end() );
497 #ifdef LIBMESH_HAVE_PETSC 508 CPPUNIT_TEST( testGhostingCouplingMatrix );
509 CPPUNIT_TEST( testGhostingNullCouplingMatrix );
510 CPPUNIT_TEST( testGhostingNullCouplingMatrixUnifRef );
512 CPPUNIT_TEST_SUITE_END();
525 this->run_ghosting_test(0,
true);
531 this->run_ghosting_test(0,
false);
537 std::unique_ptr<CouplingMatrix> coupling_matrix;
539 this->run_ghosting_test(2,
false);
546 this->build_quad_mesh(n_refinements);
549 std::unique_ptr<CouplingMatrix> coupling_matrix;
550 if (build_coupling_matrix)
551 this->setup_coupling_matrix(coupling_matrix);
559 functor.set_coupling_matrix(coupling_matrix);
574 std::unique_ptr<PointLocatorBase> point_locator = _mesh->sub_point_locator();
586 for (
const auto & elem : _mesh->active_local_subdomain_elements_ptr_range(2))
589 CPPUNIT_ASSERT_EQUAL(2, static_cast<int>(elem->subdomain_id()));
592 const std::vector<libMesh::Point> & qpoints = subdomain_two_context.
get_element_fe(u_var)->get_xyz();
598 std::set<subdomain_id_type> allowed_subdomains;
599 allowed_subdomains.insert(1);
605 for (
const auto & qp : qpoints )
607 const Elem * overlapping_elem = (*point_locator)( qp, &allowed_subdomains );
608 CPPUNIT_ASSERT(overlapping_elem);
611 subdomain_one_context.
pre_fe_reinit(system,overlapping_elem);
628 CPPUNIT_TEST( testSparsityCouplingMatrix );
629 CPPUNIT_TEST( testSparsityNullCouplingMatrix );
630 CPPUNIT_TEST( testSparsityNullCouplingMatrixUnifRef );
632 CPPUNIT_TEST_SUITE_END();
645 this->run_sparsity_pattern_test(0,
true);
651 this->run_sparsity_pattern_test(0,
false);
657 this->run_sparsity_pattern_test(1,
false);
664 this->build_quad_mesh(n_refinements);
667 std::unique_ptr<CouplingMatrix> coupling_matrix;
668 if (build_coupling_matrix)
669 this->setup_coupling_matrix(coupling_matrix);
677 coupling_functor.set_coupling_matrix(coupling_matrix);
703 std::unique_ptr<PointLocatorBase> point_locator = _mesh->sub_point_locator();
719 for (
const auto & elem : _mesh->active_local_subdomain_elements_ptr_range(1))
724 std::vector<dof_id_type> & rows = subdomain_one_context.
get_dof_indices();
735 for (
const auto & elem : _mesh->active_local_subdomain_elements_ptr_range(2))
738 CPPUNIT_ASSERT_EQUAL(2, static_cast<int>(elem->subdomain_id()));
740 const std::vector<libMesh::Point> & qpoints = subdomain_two_context.
get_element_fe(u_var)->get_xyz();
748 std::vector<dof_id_type> & rows = subdomain_two_context.
get_dof_indices();
757 std::set<subdomain_id_type> allowed_subdomains;
758 allowed_subdomains.insert(1);
764 for (
const auto & qp : qpoints )
766 const Elem * overlapping_elem = (*point_locator)( qp, &allowed_subdomains );
767 CPPUNIT_ASSERT(overlapping_elem);
770 subdomain_one_context.
pre_fe_reinit(system,overlapping_elem);
774 std::vector<dof_id_type> & v_indices = subdomain_one_context.
get_dof_indices(v_var);
775 std::vector<dof_id_type> columns(rows);
776 columns.insert( columns.end(), v_indices.begin(), v_indices.end() );
779 K21.
resize( rows.size(), columns.size() );
789 K12.
resize(v_indices.size(), rows.size());
802 #endif // LIBMESH_HAVE_PETSC 807 #ifdef LIBMESH_HAVE_PETSC class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
void run_coupling_functor_test(unsigned int n_refinements)
OverlappingCouplingFunctor(System &system)
void checkOverlappingPartitionerUnifRef()
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
virtual Node *& set_node(const unsigned int i)
virtual dof_id_type n_active_elem() const =0
This abstract base class defines the interface by which library code and user code can report associa...
virtual void pre_fe_reinit(const System &, const Elem *e)
Reinitializes local data vectors/matrices on the current geometric element.
void checkCouplingFunctorTriUnifRef()
void setup_coupling_matrix(std::unique_ptr< CouplingMatrix > &coupling)
The definition of the const_element_iterator struct.
void testSparsityNullCouplingMatrixUnifRef()
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and linear solvers ...
virtual std::unique_ptr< Partitioner > clone() const override
bool is_attached(SparseMatrix< Number > &matrix)
Matrices should not be attached more than once.
void testSparsityNullCouplingMatrix()
This is the base class from which all geometric element types are derived.
void checkCouplingFunctorQuad()
const Parallel::Communicator & comm() const
Order default_quadrature_order() const
std::unique_ptr< CouplingMatrix > _coupling_matrix
std::map< const Elem *, const CouplingMatrix *, CompareDofObjectsByPIDAndThenID > map_type
What elements do we care about and what variables do we care about on each element?
void checkCouplingFunctorQuadUnifRef()
The libMesh namespace provides an interface to certain functionality in the library.
dof_id_type n_local_dofs() const
void run_partitioner_test(unsigned int n_refinements)
const MeshBase & get_mesh() const
Real distance(const Point &p)
void add_coupling_functor(GhostingFunctor &coupling_functor, bool to_mesh=true)
Adds a functor which can specify coupling requirements for creation of sparse matrices.
void testGhostingCouplingMatrix()
dof_id_type n_dofs() const
This is the MeshBase class.
virtual void elem_fe_reinit(const std::vector< Point > *const pts=nullptr)
Reinitializes interior FE objects on the current geometric element.
The Partitioner class provides a uniform interface for partitioning algorithms.
unsigned int variable_number(std::string_view var) const
std::unique_ptr< MeshBase > _mesh
std::unique_ptr< EquationSystems > _es
virtual void _do_partition(MeshBase &mesh, const unsigned int n) override
Partition the MeshBase into n subdomains.
virtual ~OverlappingCouplingFunctor()
std::unique_ptr< PointLocatorBase > _point_locator
void testGhostingNullCouplingMatrix()
This class handles the numbering of degrees of freedom on a mesh.
void build_quad_mesh(unsigned int n_refinements=0)
void checkCouplingFunctorTri()
uint8_t processor_id_type
void reinit_send_list(MeshBase &mesh)
Clears the _send_list vector and then rebuilds it.
processor_id_type n_processors() const
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
Add the full matrix dm to the SparseMatrix.
Implements (adaptive) mesh refinement algorithms for a MeshBase.
CPPUNIT_TEST_SUITE_REGISTRATION(OverlappingFunctorTest)
static const processor_id_type invalid_processor_id
An invalid processor_id to distinguish DoFs that have not been assigned to a processor.
void assign_proc_id_subdomain(MeshBase &mesh, const subdomain_id_type sid, const dof_id_type blksize, const unsigned int n_parts, const processor_id_type offset)
Manages consistently variables, degrees of freedom, and coefficient vectors.
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
Helper function that allows us to treat a homogenous pair as a range.
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
This class provides all data required for a physics package (e.g.
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
unsigned int add_variable(std::string_view var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=nullptr)
Adds the variable var to the list of variables for this system.
void run_ghosting_test(const unsigned int n_refinements, bool build_coupling_matrix)
std::vector< T > & get_values()
void testGhostingNullCouplingMatrixUnifRef()
virtual void close()=0
Calls the SparseMatrix's internal assembly routines, ensuring that the values are consistent across p...
void set_coupling_matrix(std::unique_ptr< CouplingMatrix > &coupling_matrix)
subdomain_id_type subdomain_id() const
static std::unique_ptr< Elem > build_with_id(const ElemType type, dof_id_type id)
Calls the build() method above with a nullptr parent, and additionally sets the newly-created Elem's ...
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
void checkOverlappingPartitioner()
static std::unique_ptr< NumericVector< T > > build(const Parallel::Communicator &comm, SolverPackage solver_package=libMesh::default_solver_package(), ParallelType parallel_type=AUTOMATIC)
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
void clear_sparsity()
Clears the sparsity pattern.
This class implements specific orders of Gauss quadrature.
unsigned int level ElemType type std::set< subdomain_id_type > ss processor_id_type pid unsigned int level std::set< subdomain_id_type > virtual ss SimpleRange< element_iterator > active_subdomain_elements_ptr_range(subdomain_id_type sid)=0
std::unique_ptr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
void get_element_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for interior finite element object for variable var for the largest dimension in the mesh...
void run_sparsity_pattern_test(const unsigned int n_refinements, bool build_coupling_matrix)
unsigned int n_vars() const
void testSparsityCouplingMatrix()
void compute_sparsity(const MeshBase &)
Computes the sparsity pattern for the matrices corresponding to proc_id and sends that data to Linear...
void add_algebraic_ghosting_functor(GhostingFunctor &evaluable_functor, bool to_mesh=true)
Adds a functor which can specify algebraic ghosting requirements for use with distributed vectors...
const DofMap & get_dof_map() const
processor_id_type processor_id() const
A Point defines a location in LIBMESH_DIM dimensional Real space.
const SparseMatrix< Number > & get_system_matrix() const
const std::vector< dof_id_type > & get_send_list() const
virtual void init(const numeric_index_type m, const numeric_index_type n, const numeric_index_type m_l, const numeric_index_type n_l, const numeric_index_type nnz=30, const numeric_index_type noz=10, const numeric_index_type blocksize=1)=0
Initialize SparseMatrix with the specified sizes.
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.