From b175c31d17604aae3f807a805aee68c05b36eeb2 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sat, 19 Nov 2022 11:14:05 +0000 Subject: [PATCH 01/11] Template Geometry class over type. --- cpp/dolfinx/fem/assemble_matrix_impl.h | 17 +++++++----- cpp/dolfinx/fem/assemble_scalar_impl.h | 27 +++++++++---------- cpp/dolfinx/fem/assemble_vector_impl.h | 37 +++++++++++++------------- cpp/dolfinx/geometry/utils.cpp | 4 +-- cpp/dolfinx/io/ADIOS2Writers.cpp | 8 +++--- cpp/dolfinx/io/VTKFile.cpp | 4 +-- cpp/dolfinx/io/XDMFFile.cpp | 2 +- cpp/dolfinx/io/XDMFFile.h | 3 ++- cpp/dolfinx/io/xdmf_mesh.cpp | 4 +-- cpp/dolfinx/io/xdmf_mesh.h | 5 ++-- cpp/dolfinx/mesh/Geometry.cpp | 31 +++------------------ cpp/dolfinx/mesh/Geometry.h | 27 ++++++++++++------- cpp/dolfinx/mesh/Mesh.cpp | 8 +++--- cpp/dolfinx/mesh/Mesh.h | 11 ++++---- cpp/dolfinx/mesh/utils.cpp | 2 +- python/dolfinx/wrappers/mesh.cpp | 20 +++++++------- 16 files changed, 100 insertions(+), 110 deletions(-) diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h index e5aefd0bf08..9262d8068dd 100644 --- a/cpp/dolfinx/fem/assemble_matrix_impl.h +++ b/cpp/dolfinx/fem/assemble_matrix_impl.h @@ -27,7 +27,8 @@ namespace dolfinx::fem::impl /// Execute kernel over cells and accumulate result in matrix template void assemble_cells( - la::MatSet auto mat_set, const mesh::Geometry& geometry, + la::MatSet auto mat_set, + const mesh::Geometry> geometry, std::span cells, const std::function&, const std::span&, @@ -47,7 +48,7 @@ void assemble_cells( // Prepare cell geometry const graph::AdjacencyList& x_dofmap = geometry.dofmap(); const std::size_t num_dofs_g = geometry.cmap().dim(); - std::span x = geometry.x(); + auto x = geometry.x(); // Iterate over active cells const int num_dofs0 = dofmap0.links(0).size(); @@ -122,7 +123,8 @@ void assemble_cells( /// Execute kernel over exterior facets and accumulate result in Mat template void assemble_exterior_facets( - la::MatSet auto mat_set, const mesh::Geometry& geometry, + la::MatSet auto mat_set, + const mesh::Geometry>& geometry, std::span facets, const std::function&, const std::span&, @@ -142,7 +144,7 @@ void assemble_exterior_facets( // Prepare cell geometry const graph::AdjacencyList& x_dofmap = geometry.dofmap(); const std::size_t num_dofs_g = geometry.cmap().dim(); - std::span x = geometry.x(); + auto x = geometry.x(); // Data structures used in assembly std::vector> coordinate_dofs(3 * num_dofs_g); @@ -216,8 +218,9 @@ void assemble_exterior_facets( /// Execute kernel over interior facets and accumulate result in Mat template void assemble_interior_facets( - la::MatSet auto mat_set, const mesh::Geometry& geometry, - int num_cell_facets, std::span facets, + la::MatSet auto mat_set, + const mesh::Geometry>& geometry, int num_cell_facets, + std::span facets, const std::function&, const std::span&, std::int32_t, int)>& dof_transform, @@ -237,7 +240,7 @@ void assemble_interior_facets( // Prepare cell geometry const graph::AdjacencyList& x_dofmap = geometry.dofmap(); const std::size_t num_dofs_g = geometry.cmap().dim(); - std::span x = geometry.x(); + auto x = geometry.x(); // Data structures used in assembly using X = scalar_value_type_t; diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index a968f0986ac..a6ab68b6941 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -23,7 +23,7 @@ namespace dolfinx::fem::impl /// Assemble functional over cells template -T assemble_cells(const mesh::Geometry& geometry, +T assemble_cells(const mesh::Geometry>& geometry, std::span cells, FEkernel auto fn, std::span constants, std::span coeffs, int cstride) @@ -35,7 +35,7 @@ T assemble_cells(const mesh::Geometry& geometry, // Prepare cell geometry const graph::AdjacencyList& x_dofmap = geometry.dofmap(); const std::size_t num_dofs_g = geometry.cmap().dim(); - std::span x = geometry.x(); + auto x = geometry.x(); // Create data structures used in assembly std::vector> coordinate_dofs(3 * num_dofs_g); @@ -63,10 +63,10 @@ T assemble_cells(const mesh::Geometry& geometry, /// Execute kernel over exterior facets and accumulate result template -T assemble_exterior_facets(const mesh::Geometry& geometry, - std::span facets, - FEkernel auto fn, std::span constants, - std::span coeffs, int cstride) +T assemble_exterior_facets( + const mesh::Geometry>& geometry, + std::span facets, FEkernel auto fn, + std::span constants, std::span coeffs, int cstride) { T value(0); if (facets.empty()) @@ -75,7 +75,7 @@ T assemble_exterior_facets(const mesh::Geometry& geometry, // Prepare cell geometry const graph::AdjacencyList& x_dofmap = geometry.dofmap(); const std::size_t num_dofs_g = geometry.cmap().dim(); - std::span x = geometry.x(); + auto x = geometry.x(); // Create data structures used in assembly std::vector> coordinate_dofs(3 * num_dofs_g); @@ -105,12 +105,11 @@ T assemble_exterior_facets(const mesh::Geometry& geometry, /// Assemble functional over interior facets template -T assemble_interior_facets(const mesh::Geometry& geometry, int num_cell_facets, - std::span facets, - FEkernel auto fn, std::span constants, - std::span coeffs, int cstride, - std::span offsets, - std::span perms) +T assemble_interior_facets( + const mesh::Geometry>& geometry, int num_cell_facets, + std::span facets, FEkernel auto fn, + std::span constants, std::span coeffs, int cstride, + std::span offsets, std::span perms) { T value(0); if (facets.empty()) @@ -119,7 +118,7 @@ T assemble_interior_facets(const mesh::Geometry& geometry, int num_cell_facets, // Prepare cell geometry const graph::AdjacencyList& x_dofmap = geometry.dofmap(); const std::size_t num_dofs_g = geometry.cmap().dim(); - std::span x = geometry.x(); + auto x = geometry.x(); // Create data structures used in assembly using X = scalar_value_type_t; diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index 65cb679e931..bb96d9aeea4 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -37,8 +37,8 @@ namespace dolfinx::fem::impl /// @tparam _bs1 The block size of the trial function dof map. template void _lift_bc_cells( - std::span b, const mesh::Geometry& geometry, FEkernel auto kernel, - std::span cells, + std::span b, const mesh::Geometry>& geometry, + FEkernel auto kernel, std::span cells, const std::function&, const std::span&, std::int32_t, int)>& dof_transform, @@ -61,7 +61,7 @@ void _lift_bc_cells( // Prepare cell geometry const graph::AdjacencyList& x_dofmap = geometry.dofmap(); const std::size_t num_dofs_g = geometry.cmap().dim(); - std::span x = geometry.x(); + auto x = geometry.x(); // Data structures used in bc application std::vector> coordinate_dofs(3 * num_dofs_g); @@ -194,8 +194,8 @@ void _lift_bc_cells( /// @tparam _bs1 The block size of the trial function dof map. template void _lift_bc_exterior_facets( - std::span b, const mesh::Geometry& geometry, FEkernel auto kernel, - std::span facets, + std::span b, const mesh::Geometry>& geometry, + FEkernel auto kernel, std::span facets, const std::function&, const std::span&, std::int32_t, int)>& dof_transform, @@ -215,7 +215,7 @@ void _lift_bc_exterior_facets( // Prepare cell geometry const graph::AdjacencyList& x_dofmap = geometry.dofmap(); const std::size_t num_dofs_g = geometry.cmap().dim(); - std::span x = geometry.x(); + auto x = geometry.x(); // Data structures used in bc application std::vector> coordinate_dofs(3 * num_dofs_g); @@ -302,8 +302,9 @@ void _lift_bc_exterior_facets( /// @tparam _bs1 The block size of the trial function dof map. template void _lift_bc_interior_facets( - std::span b, const mesh::Geometry& geometry, int num_cell_facets, - FEkernel auto kernel, std::span facets, + std::span b, const mesh::Geometry>& geometry, + int num_cell_facets, FEkernel auto kernel, + std::span facets, const std::function&, const std::span&, std::int32_t, int)>& dof_transform, @@ -324,7 +325,7 @@ void _lift_bc_interior_facets( // Prepare cell geometry const graph::AdjacencyList& x_dofmap = geometry.dofmap(); const std::size_t num_dofs_g = geometry.cmap().dim(); - std::span x = geometry.x(); + auto x = geometry.x(); // Data structures used in assembly using X = scalar_value_type_t; @@ -490,7 +491,7 @@ void assemble_cells( const std::function&, const std::span&, std::int32_t, int)>& dof_transform, - std::span b, const mesh::Geometry& geometry, + std::span b, const mesh::Geometry>& geometry, std::span cells, const graph::AdjacencyList& dofmap, int bs, FEkernel auto kernel, std::span constants, @@ -505,7 +506,7 @@ void assemble_cells( // Prepare cell geometry const graph::AdjacencyList& x_dofmap = geometry.dofmap(); const std::size_t num_dofs_g = geometry.cmap().dim(); - std::span x = geometry.x(); + auto x = geometry.x(); // FIXME: Add proper interface for num_dofs // Create data structures used in assembly @@ -561,7 +562,7 @@ void assemble_exterior_facets( const std::function&, const std::span&, std::int32_t, int)>& dof_transform, - std::span b, const mesh::Geometry& geometry, + std::span b, const mesh::Geometry>& geometry, std::span facets, const graph::AdjacencyList& dofmap, int bs, FEkernel auto fn, std::span constants, @@ -576,7 +577,7 @@ void assemble_exterior_facets( // Prepare cell geometry const graph::AdjacencyList& x_dofmap = geometry.dofmap(); const std::size_t num_dofs_g = geometry.cmap().dim(); - std::span x = geometry.x(); + auto x = geometry.x(); // FIXME: Add proper interface for num_dofs // Create data structures used in assembly @@ -633,17 +634,17 @@ void assemble_interior_facets( const std::function&, const std::span&, std::int32_t, int)>& dof_transform, - std::span b, const mesh::Geometry& geometry, int num_cell_facets, - std::span facets, const fem::DofMap& dofmap, - FEkernel auto fn, std::span constants, - std::span coeffs, int cstride, + std::span b, const mesh::Geometry>& geometry, + int num_cell_facets, std::span facets, + const fem::DofMap& dofmap, FEkernel auto fn, + std::span constants, std::span coeffs, int cstride, std::span cell_info, const std::function& get_perm) { // Prepare cell geometry const graph::AdjacencyList& x_dofmap = geometry.dofmap(); const std::size_t num_dofs_g = geometry.cmap().dim(); - std::span x = geometry.x(); + auto x = geometry.x(); // Create data structures used in assembly using X = scalar_value_type_t; diff --git a/cpp/dolfinx/geometry/utils.cpp b/cpp/dolfinx/geometry/utils.cpp index 691bc672bea..cdacedc4724 100644 --- a/cpp/dolfinx/geometry/utils.cpp +++ b/cpp/dolfinx/geometry/utils.cpp @@ -392,7 +392,7 @@ geometry::shortest_vector(const mesh::Mesh& mesh, int dim, std::span points) { const int tdim = mesh.topology().dim(); - const mesh::Geometry& geometry = mesh.geometry(); + const mesh::Geometry& geometry = mesh.geometry(); std::span geom_dofs = geometry.x(); const graph::AdjacencyList& x_dofmap = geometry.dofmap(); std::vector shortest_vectors(3 * entities.size()); @@ -515,7 +515,7 @@ int geometry::compute_first_colliding_cell( else { constexpr double eps2 = 1e-20; - const mesh::Geometry& geometry = mesh.geometry(); + const mesh::Geometry& geometry = mesh.geometry(); std::span geom_dofs = geometry.x(); const graph::AdjacencyList& x_dofmap = geometry.dofmap(); const std::size_t num_nodes = geometry.cmap().dim(); diff --git a/cpp/dolfinx/io/ADIOS2Writers.cpp b/cpp/dolfinx/io/ADIOS2Writers.cpp index ec231fc4796..deea238cf27 100644 --- a/cpp/dolfinx/io/ADIOS2Writers.cpp +++ b/cpp/dolfinx/io/ADIOS2Writers.cpp @@ -145,7 +145,7 @@ void vtx_write_data(adios2::IO& io, adios2::Engine& engine, void vtx_write_mesh(adios2::IO& io, adios2::Engine& engine, const mesh::Mesh& mesh) { - const mesh::Geometry& geometry = mesh.geometry(); + const mesh::Geometry& geometry = mesh.geometry(); const mesh::Topology& topology = mesh.topology(); // "Put" geometry @@ -411,7 +411,7 @@ std::vector pack_function_data(const fem::Function& u) assert(dofmap); auto mesh = V->mesh(); assert(mesh); - const mesh::Geometry& geometry = mesh->geometry(); + const mesh::Geometry& geometry = mesh->geometry(); const mesh::Topology& topology = mesh->topology(); // The Function and the mesh must have identical element_dof_layouts @@ -534,7 +534,7 @@ void fides_write_data(adios2::IO& io, adios2::Engine& engine, void fides_write_mesh(adios2::IO& io, adios2::Engine& engine, const mesh::Mesh& mesh) { - const mesh::Geometry& geometry = mesh.geometry(); + const mesh::Geometry& geometry = mesh.geometry(); const mesh::Topology& topology = mesh.topology(); // "Put" geometry data @@ -568,7 +568,7 @@ void fides_write_mesh(adios2::IO& io, adios2::Engine& engine, /// @param[in] mesh The mesh void fides_initialize_mesh_attributes(adios2::IO& io, const mesh::Mesh& mesh) { - const mesh::Geometry& geometry = mesh.geometry(); + const mesh::Geometry& geometry = mesh.geometry(); const mesh::Topology& topology = mesh.topology(); // Check that mesh is first order mesh diff --git a/cpp/dolfinx/io/VTKFile.cpp b/cpp/dolfinx/io/VTKFile.cpp index 413a62d8e87..fcd4036454b 100644 --- a/cpp/dolfinx/io/VTKFile.cpp +++ b/cpp/dolfinx/io/VTKFile.cpp @@ -421,7 +421,7 @@ void write_function( std::vector tmp; std::tie(tmp, cshape) = io::extract_vtk_connectivity(*mesh0); cells.assign(tmp.begin(), tmp.end()); - const mesh::Geometry& geometry = mesh0->geometry(); + const mesh::Geometry& geometry = mesh0->geometry(); x.assign(geometry.x().begin(), geometry.x().end()); xshape = {geometry.x().size() / 3, 3}; x_id = geometry.input_global_indices(); @@ -749,7 +749,7 @@ void io::VTKFile::write(const mesh::Mesh& mesh, double time) // Get mesh data for this rank const mesh::Topology& topology = mesh.topology(); - const mesh::Geometry& geometry = mesh.geometry(); + const mesh::Geometry& geometry = mesh.geometry(); auto xmap = geometry.index_map(); assert(xmap); const int tdim = topology.dim(); diff --git a/cpp/dolfinx/io/XDMFFile.cpp b/cpp/dolfinx/io/XDMFFile.cpp index d0bb5fd8292..94f8a0f884a 100644 --- a/cpp/dolfinx/io/XDMFFile.cpp +++ b/cpp/dolfinx/io/XDMFFile.cpp @@ -202,7 +202,7 @@ void XDMFFile::write_mesh(const mesh::Mesh& mesh, std::string xpath) _xml_doc->save_file(_filename.c_str(), " "); } //----------------------------------------------------------------------------- -void XDMFFile::write_geometry(const mesh::Geometry& geometry, std::string name, +void XDMFFile::write_geometry(const mesh::Geometry& geometry, std::string name, std::string xpath) { pugi::xml_node node = _xml_doc->select_node(xpath.c_str()).node(); diff --git a/cpp/dolfinx/io/XDMFFile.h b/cpp/dolfinx/io/XDMFFile.h index b53cf95c735..182964f9a63 100644 --- a/cpp/dolfinx/io/XDMFFile.h +++ b/cpp/dolfinx/io/XDMFFile.h @@ -32,6 +32,7 @@ class Function; namespace dolfinx::mesh { +template class Geometry; enum class GhostMode : int; class Mesh; @@ -89,7 +90,7 @@ class XDMFFile /// @param[in] geometry /// @param[in] name /// @param[in] xpath XPath of a node where Geometry will be inserted - void write_geometry(const mesh::Geometry& geometry, std::string name, + void write_geometry(const mesh::Geometry& geometry, std::string name, std::string xpath = "/Xdmf/Domain"); /// Read in Mesh diff --git a/cpp/dolfinx/io/xdmf_mesh.cpp b/cpp/dolfinx/io/xdmf_mesh.cpp index 66663ab611a..2f9cf861491 100644 --- a/cpp/dolfinx/io/xdmf_mesh.cpp +++ b/cpp/dolfinx/io/xdmf_mesh.cpp @@ -21,7 +21,7 @@ void xdmf_mesh::add_topology_data(MPI_Comm comm, pugi::xml_node& xml_node, const hid_t h5_id, const std::string path_prefix, const mesh::Topology& topology, - const mesh::Geometry& geometry, int dim, + const mesh::Geometry& geometry, int dim, std::span entities) { LOG(INFO) << "Adding topology data to node \"" << xml_node.path('/') << "\""; @@ -149,7 +149,7 @@ void xdmf_mesh::add_topology_data(MPI_Comm comm, pugi::xml_node& xml_node, void xdmf_mesh::add_geometry_data(MPI_Comm comm, pugi::xml_node& xml_node, const hid_t h5_id, const std::string path_prefix, - const mesh::Geometry& geometry) + const mesh::Geometry& geometry) { LOG(INFO) << "Adding geometry data to node \"" << xml_node.path('/') << "\""; auto map = geometry.index_map(); diff --git a/cpp/dolfinx/io/xdmf_mesh.h b/cpp/dolfinx/io/xdmf_mesh.h index 2e4f150013b..1c6a4391168 100644 --- a/cpp/dolfinx/io/xdmf_mesh.h +++ b/cpp/dolfinx/io/xdmf_mesh.h @@ -24,6 +24,7 @@ namespace dolfinx namespace mesh { +template class Geometry; class Mesh; class Topology; @@ -53,13 +54,13 @@ void add_mesh(MPI_Comm comm, pugi::xml_node& xml_node, const hid_t h5_id, void add_topology_data(MPI_Comm comm, pugi::xml_node& xml_node, const hid_t h5_id, const std::string path_prefix, const mesh::Topology& topology, - const mesh::Geometry& geometry, int cell_dim, + const mesh::Geometry& geometry, int cell_dim, std::span entities); /// Add Geometry xml node void add_geometry_data(MPI_Comm comm, pugi::xml_node& xml_node, const hid_t h5_id, const std::string path_prefix, - const mesh::Geometry& geometry); + const mesh::Geometry& geometry); /// @brief Read geometry (coordinate) data. /// diff --git a/cpp/dolfinx/mesh/Geometry.cpp b/cpp/dolfinx/mesh/Geometry.cpp index ffc19cf02a4..cf2d1038c09 100644 --- a/cpp/dolfinx/mesh/Geometry.cpp +++ b/cpp/dolfinx/mesh/Geometry.cpp @@ -16,32 +16,7 @@ using namespace dolfinx; using namespace dolfinx::mesh; //----------------------------------------------------------------------------- -int Geometry::dim() const { return _dim; } -//----------------------------------------------------------------------------- -const graph::AdjacencyList& Geometry::dofmap() const -{ - return _dofmap; -} -//----------------------------------------------------------------------------- -std::shared_ptr Geometry::index_map() const -{ - return _index_map; -} -//----------------------------------------------------------------------------- -std::span Geometry::x() { return _x; } -//----------------------------------------------------------------------------- -std::span Geometry::x() const { return _x; } -//----------------------------------------------------------------------------- -const fem::CoordinateElement& Geometry::cmap() const { return _cmap; } -//----------------------------------------------------------------------------- -const std::vector& Geometry::input_global_indices() const -{ - return _input_global_indices; -} -//----------------------------------------------------------------------------- - -//----------------------------------------------------------------------------- -mesh::Geometry mesh::create_geometry( +mesh::Geometry mesh::create_geometry( MPI_Comm comm, const Topology& topology, const fem::CoordinateElement& element, const graph::AdjacencyList& cell_nodes, @@ -116,7 +91,7 @@ mesh::Geometry mesh::create_geometry( std::next(xg.begin(), 3 * i)); } - return Geometry(dof_index_map, std::move(dofmap), element, std::move(xg), dim, - std::move(igi)); + return Geometry(dof_index_map, std::move(dofmap), element, + std::move(xg), dim, std::move(igi)); } //----------------------------------------------------------------------------- diff --git a/cpp/dolfinx/mesh/Geometry.h b/cpp/dolfinx/mesh/Geometry.h index 9dac3a37b4b..d50d5155598 100644 --- a/cpp/dolfinx/mesh/Geometry.h +++ b/cpp/dolfinx/mesh/Geometry.h @@ -25,6 +25,7 @@ namespace dolfinx::mesh class Topology; /// @brief Geometry stores the geometry imposed on a mesh. +template class Geometry { public: @@ -41,7 +42,7 @@ class Geometry /// point, commonly from a mesh input file. The type is /// `std:vector`. template > U, - std::convertible_to> V, + std::convertible_to> V, std::convertible_to> W> Geometry(std::shared_ptr index_map, U&& dofmap, const fem::CoordinateElement& element, V&& x, int dim, @@ -71,34 +72,40 @@ class Geometry Geometry& operator=(Geometry&&) = default; /// Return Euclidean dimension of coordinate system - int dim() const; + int dim() const { return _dim; } /// DOF map - const graph::AdjacencyList& dofmap() const; + const graph::AdjacencyList& dofmap() const { return _dofmap; } /// Index map - std::shared_ptr index_map() const; + std::shared_ptr index_map() const + { + return _index_map; + } /// @brief Access geometry degrees-of-freedom data (const version). /// /// @return The flattened row-major geometry data, where the shape is /// (num_points, 3) - std::span x() const; + std::span x() const { return _x; } /// @brief Access geometry degrees-of-freedom data (non-const /// version). /// /// @return The flattened row-major geometry data, where the shape is /// (num_points, 3) - std::span x(); + std::span x() { return _x; } /// @brief The element that describes the geometry map. /// /// @return The coordinate/geometry element - const fem::CoordinateElement& cmap() const; + const fem::CoordinateElement& cmap() const { return _cmap; } /// Global user indices - const std::vector& input_global_indices() const; + const std::vector& input_global_indices() const + { + return _input_global_indices; + } private: // Geometric dimension @@ -115,7 +122,7 @@ class Geometry // Coordinates for all points stored as a contiguous array (row-major, // column size = 3) - std::vector _x; + std::vector _x; // Global indices as provided on Geometry creation std::vector _input_global_indices; @@ -141,7 +148,7 @@ class Geometry /// @param[in] dim The geometric dimension (1, 2, or 3) /// @param[in] reorder_fn Function for re-ordering the degree-of-freedom /// map associated with the geometry data -mesh::Geometry +mesh::Geometry create_geometry(MPI_Comm comm, const Topology& topology, const fem::CoordinateElement& element, const graph::AdjacencyList& cells, diff --git a/cpp/dolfinx/mesh/Mesh.cpp b/cpp/dolfinx/mesh/Mesh.cpp index 277b56c3eb5..eaf89e4db8c 100644 --- a/cpp/dolfinx/mesh/Mesh.cpp +++ b/cpp/dolfinx/mesh/Mesh.cpp @@ -314,7 +314,7 @@ mesh::create_submesh(const Mesh& mesh, int dim, submesh_topology.set_connectivity(submesh_e_to_v, dim, 0); // -- Submesh geometry - const Geometry& geometry = mesh.geometry(); + const Geometry& geometry = mesh.geometry(); // Get the geometry dofs in the submesh based on the entities in // submesh @@ -454,7 +454,7 @@ mesh::create_submesh(const Mesh& mesh, int dim, { return mesh_igi[submesh_x_dof]; }); // Create geometry - Geometry submesh_geometry( + Geometry submesh_geometry( submesh_x_dof_index_map, std::move(submesh_x_dofmap), submesh_coord_ele, std::move(submesh_x), geometry.dim(), std::move(submesh_igi)); @@ -470,9 +470,9 @@ const Topology& Mesh::topology() const { return _topology; } //----------------------------------------------------------------------------- Topology& Mesh::topology_mutable() const { return _topology; } //----------------------------------------------------------------------------- -Geometry& Mesh::geometry() { return _geometry; } +Geometry& Mesh::geometry() { return _geometry; } //----------------------------------------------------------------------------- -const Geometry& Mesh::geometry() const { return _geometry; } +const Geometry& Mesh::geometry() const { return _geometry; } //----------------------------------------------------------------------------- MPI_Comm Mesh::comm() const { return _comm.comm(); } //----------------------------------------------------------------------------- diff --git a/cpp/dolfinx/mesh/Mesh.h b/cpp/dolfinx/mesh/Mesh.h index 249f121dc67..2305bba3f96 100644 --- a/cpp/dolfinx/mesh/Mesh.h +++ b/cpp/dolfinx/mesh/Mesh.h @@ -37,7 +37,8 @@ class Mesh /// @param[in] comm MPI Communicator /// @param[in] topology Mesh topology /// @param[in] geometry Mesh geometry - template U, std::convertible_to V> + template U, + std::convertible_to> V> Mesh(MPI_Comm comm, U&& topology, V&& geometry) : _topology(std::forward(topology)), _geometry(std::forward(geometry)), _comm(comm) @@ -80,17 +81,17 @@ class Mesh /// Get mesh geometry /// @return The geometry object associated with the mesh - Geometry& geometry(); + Geometry& geometry(); /// Get mesh geometry (const version) /// @return The geometry object associated with the mesh - const Geometry& geometry() const; + const Geometry& geometry() const; /// Mesh MPI communicator /// @return The communicator on which the mesh is distributed MPI_Comm comm() const; - /// Name + /// Nameg std::string name = "mesh"; private: @@ -102,7 +103,7 @@ class Mesh mutable Topology _topology; // Mesh geometry - Geometry _geometry; + Geometry _geometry; // MPI communicator dolfinx::MPI::Comm _comm; diff --git a/cpp/dolfinx/mesh/utils.cpp b/cpp/dolfinx/mesh/utils.cpp index 6e63fa4a86f..c9e1bbfcdf3 100644 --- a/cpp/dolfinx/mesh/utils.cpp +++ b/cpp/dolfinx/mesh/utils.cpp @@ -500,7 +500,7 @@ mesh::entities_to_geometry(const Mesh& mesh, int dim, if (orient and (cell_type != CellType::tetrahedron or dim != 2)) throw std::runtime_error("Can only orient facets of a tetrahedral mesh"); - const Geometry& geometry = mesh.geometry(); + const Geometry& geometry = mesh.geometry(); auto x = geometry.x(); const Topology& topology = mesh.topology(); diff --git a/python/dolfinx/wrappers/mesh.cpp b/python/dolfinx/wrappers/mesh.cpp index 56f9e2a5f46..2370eab7755 100644 --- a/python/dolfinx/wrappers/mesh.cpp +++ b/python/dolfinx/wrappers/mesh.cpp @@ -250,25 +250,27 @@ void mesh(py::module& m) .value("shared_vertex", dolfinx::mesh::GhostMode::shared_vertex); // dolfinx::mesh::Geometry class - py::class_>( + py::class_, + std::shared_ptr>>( m, "Geometry", "Geometry object") - .def_property_readonly("dim", &dolfinx::mesh::Geometry::dim, + .def_property_readonly("dim", &dolfinx::mesh::Geometry::dim, "Geometric dimension") - .def_property_readonly("dofmap", &dolfinx::mesh::Geometry::dofmap) - .def("index_map", &dolfinx::mesh::Geometry::index_map) + .def_property_readonly("dofmap", &dolfinx::mesh::Geometry::dofmap) + .def("index_map", &dolfinx::mesh::Geometry::index_map) .def_property_readonly( "x", - [](const dolfinx::mesh::Geometry& self) + [](const dolfinx::mesh::Geometry& self) { std::array shape = {self.x().size() / 3, 3}; return py::array_t(shape, self.x().data(), py::cast(self)); }, "Return coordinates of all geometry points. Each row is the " "coordinate of a point.") - .def_property_readonly("cmap", &dolfinx::mesh::Geometry::cmap, + .def_property_readonly("cmap", &dolfinx::mesh::Geometry::cmap, "The coordinate map") - .def_property_readonly("input_global_indices", - &dolfinx::mesh::Geometry::input_global_indices); + .def_property_readonly( + "input_global_indices", + &dolfinx::mesh::Geometry::input_global_indices); // dolfinx::mesh::TopologyComputation m.def( @@ -340,7 +342,7 @@ void mesh(py::module& m) .def(py::init( [](const MPICommWrapper comm, const dolfinx::mesh::Topology& topology, - dolfinx::mesh::Geometry& geometry) + dolfinx::mesh::Geometry& geometry) { return dolfinx::mesh::Mesh(comm.get(), topology, geometry); }), py::arg("comm"), py::arg("topology"), py::arg("geometry")) .def_property_readonly( From 451ecec701cdd4912a49d2516ec6976656b9c9c1 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sat, 19 Nov 2022 19:20:39 +0000 Subject: [PATCH 02/11] Work on Geometry type --- cpp/dolfinx/fem/assemble_scalar_impl.h | 51 ++++++++++++----- cpp/dolfinx/fem/assemble_vector_impl.h | 77 +++++++++++++++++--------- cpp/dolfinx/mesh/Geometry.h | 19 +++++-- 3 files changed, 101 insertions(+), 46 deletions(-) diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index a6ab68b6941..21a2b3bbec8 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -23,10 +23,10 @@ namespace dolfinx::fem::impl /// Assemble functional over cells template -T assemble_cells(const mesh::Geometry>& geometry, - std::span cells, FEkernel auto fn, - std::span constants, std::span coeffs, - int cstride) +T assemble_cells( + const mesh::Geometry>>& geometry, + std::span cells, FEkernel auto fn, + std::span constants, std::span coeffs, int cstride) { T value(0); if (cells.empty()) @@ -64,7 +64,7 @@ T assemble_cells(const mesh::Geometry>& geometry, /// Execute kernel over exterior facets and accumulate result template T assemble_exterior_facets( - const mesh::Geometry>& geometry, + const mesh::Geometry>& geometry, std::span facets, FEkernel auto fn, std::span constants, std::span coeffs, int cstride) { @@ -106,7 +106,7 @@ T assemble_exterior_facets( /// Assemble functional over interior facets template T assemble_interior_facets( - const mesh::Geometry>& geometry, int num_cell_facets, + const mesh::Geometry>& geometry, int num_cell_facets, std::span facets, FEkernel auto fn, std::span constants, std::span coeffs, int cstride, std::span offsets, std::span perms) @@ -160,10 +160,12 @@ T assemble_interior_facets( return value; } -/// Assemble functional into an scalar +/// Assemble functional into an scalar with provided mesh geometry. template T assemble_scalar( - const fem::Form& M, std::span constants, + const fem::Form& M, + const mesh::Geometry>& geometry, + std::span constants, const std::map, std::pair, int>>& coefficients) { @@ -176,8 +178,8 @@ T assemble_scalar( const auto& fn = M.kernel(IntegralType::cell, i); const auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i}); const std::vector& cells = M.cell_domains(i); - value += impl::assemble_cells(mesh->geometry(), cells, fn, constants, - coeffs, cstride); + value += impl::assemble_cells(geometry, cells, fn, constants, coeffs, + cstride); } for (int i : M.integral_ids(IntegralType::exterior_facet)) @@ -186,8 +188,8 @@ T assemble_scalar( const auto& [coeffs, cstride] = coefficients.at({IntegralType::exterior_facet, i}); const std::vector& facets = M.exterior_facet_domains(i); - value += impl::assemble_exterior_facets(mesh->geometry(), facets, fn, - constants, coeffs, cstride); + value += impl::assemble_exterior_facets(geometry, facets, fn, constants, + coeffs, cstride); } if (M.num_integrals(IntegralType::interior_facet) > 0) @@ -206,13 +208,32 @@ T assemble_scalar( const auto& [coeffs, cstride] = coefficients.at({IntegralType::interior_facet, i}); const std::vector& facets = M.interior_facet_domains(i); - value += impl::assemble_interior_facets(mesh->geometry(), num_cell_facets, - facets, fn, constants, coeffs, - cstride, c_offsets, perms); + value += impl::assemble_interior_facets(geometry, num_cell_facets, facets, + fn, constants, coeffs, cstride, + c_offsets, perms); } } return value; } +/// Assemble functional into an scalar +template +T assemble_scalar( + const fem::Form& M, std::span constants, + const std::map, + std::pair, int>>& coefficients) +{ + std::shared_ptr mesh = M.mesh(); + assert(mesh); + if constexpr (std::is_same_v>) + return assemble_scalar(M, mesh->geometry(), constants, coefficients); + else + { + + return assemble_scalar(M, mesh->geometry().astype>(), + constants, coefficients); + } +} + } // namespace dolfinx::fem::impl diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index bb96d9aeea4..83e79ee3e9a 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -491,7 +491,7 @@ void assemble_cells( const std::function&, const std::span&, std::int32_t, int)>& dof_transform, - std::span b, const mesh::Geometry>& geometry, + std::span b, const mesh::Geometry>& geometry, std::span cells, const graph::AdjacencyList& dofmap, int bs, FEkernel auto kernel, std::span constants, @@ -562,7 +562,7 @@ void assemble_exterior_facets( const std::function&, const std::span&, std::int32_t, int)>& dof_transform, - std::span b, const mesh::Geometry>& geometry, + std::span b, const mesh::Geometry>& geometry, std::span facets, const graph::AdjacencyList& dofmap, int bs, FEkernel auto fn, std::span constants, @@ -634,7 +634,7 @@ void assemble_interior_facets( const std::function&, const std::span&, std::int32_t, int)>& dof_transform, - std::span b, const mesh::Geometry>& geometry, + std::span b, const mesh::Geometry>& geometry, int num_cell_facets, std::span facets, const fem::DofMap& dofmap, FEkernel auto fn, std::span constants, std::span coeffs, int cstride, @@ -934,7 +934,9 @@ void apply_lifting( /// @param[in] coefficients Packed coefficients that appear in `L` template void assemble_vector( - std::span b, const Form& L, std::span constants, + std::span b, const Form& L, + const mesh::Geometry>& geometry, + std::span constants, const std::map, std::pair, int>>& coefficients) { @@ -972,20 +974,18 @@ void assemble_vector( const std::vector& cells = L.cell_domains(i); if (bs == 1) { - impl::assemble_cells(dof_transform, b, mesh->geometry(), cells, - dofs, bs, fn, constants, coeffs, cstride, - cell_info); + impl::assemble_cells(dof_transform, b, geometry, cells, dofs, bs, + fn, constants, coeffs, cstride, cell_info); } else if (bs == 3) { - impl::assemble_cells(dof_transform, b, mesh->geometry(), cells, - dofs, bs, fn, constants, coeffs, cstride, - cell_info); + impl::assemble_cells(dof_transform, b, geometry, cells, dofs, bs, + fn, constants, coeffs, cstride, cell_info); } else { - impl::assemble_cells(dof_transform, b, mesh->geometry(), cells, dofs, bs, - fn, constants, coeffs, cstride, cell_info); + impl::assemble_cells(dof_transform, b, geometry, cells, dofs, bs, fn, + constants, coeffs, cstride, cell_info); } } @@ -997,20 +997,20 @@ void assemble_vector( const std::vector& facets = L.exterior_facet_domains(i); if (bs == 1) { - impl::assemble_exterior_facets(dof_transform, b, mesh->geometry(), - facets, dofs, bs, fn, constants, - coeffs, cstride, cell_info); + impl::assemble_exterior_facets(dof_transform, b, geometry, facets, + dofs, bs, fn, constants, coeffs, + cstride, cell_info); } else if (bs == 3) { - impl::assemble_exterior_facets(dof_transform, b, mesh->geometry(), - facets, dofs, bs, fn, constants, - coeffs, cstride, cell_info); + impl::assemble_exterior_facets(dof_transform, b, geometry, facets, + dofs, bs, fn, constants, coeffs, + cstride, cell_info); } else { - impl::assemble_exterior_facets(dof_transform, b, mesh->geometry(), facets, - dofs, bs, fn, constants, coeffs, cstride, + impl::assemble_exterior_facets(dof_transform, b, geometry, facets, dofs, + bs, fn, constants, coeffs, cstride, cell_info); } } @@ -1039,22 +1039,47 @@ void assemble_vector( if (bs == 1) { impl::assemble_interior_facets( - dof_transform, b, mesh->geometry(), num_cell_facets, facets, - *dofmap, fn, constants, coeffs, cstride, cell_info, get_perm); + dof_transform, b, geometry, num_cell_facets, facets, *dofmap, fn, + constants, coeffs, cstride, cell_info, get_perm); } else if (bs == 3) { impl::assemble_interior_facets( - dof_transform, b, mesh->geometry(), num_cell_facets, facets, - *dofmap, fn, constants, coeffs, cstride, cell_info, get_perm); + dof_transform, b, geometry, num_cell_facets, facets, *dofmap, fn, + constants, coeffs, cstride, cell_info, get_perm); } else { impl::assemble_interior_facets( - dof_transform, b, mesh->geometry(), num_cell_facets, facets, - *dofmap, fn, constants, coeffs, cstride, cell_info, get_perm); + dof_transform, b, geometry, num_cell_facets, facets, *dofmap, fn, + constants, coeffs, cstride, cell_info, get_perm); } } } } + +/// Assemble linear form into a vector +/// @param[in,out] b The vector to be assembled. It will not be zeroed +/// before assembly. +/// @param[in] L The linear forms to assemble into b +/// @param[in] constants Packed constants that appear in `L` +/// @param[in] coefficients Packed coefficients that appear in `L` +template +void assemble_vector( + std::span b, const Form& L, std::span constants, + const std::map, + std::pair, int>>& coefficients) +{ + std::shared_ptr mesh = L.mesh(); + assert(mesh); + if constexpr (std::is_same_v>) + return assemble_vector(b, L, mesh->geometry(), constants, coefficients); + else + { + + return assemble_vector(b, L, + mesh->geometry().astype>(), + constants, coefficients); + } +} } // namespace dolfinx::fem::impl diff --git a/cpp/dolfinx/mesh/Geometry.h b/cpp/dolfinx/mesh/Geometry.h index d50d5155598..8f63aeb4814 100644 --- a/cpp/dolfinx/mesh/Geometry.h +++ b/cpp/dolfinx/mesh/Geometry.h @@ -1,4 +1,4 @@ -// Copyright (C) 2006-2020 Anders Logg and Garth N. Wells +// Copyright (C) 2006-2022 Anders Logg and Garth N. Wells // // This file is part of DOLFINx (https://www.fenicsproject.org) // @@ -35,8 +35,8 @@ class Geometry /// @param[in] dofmap The geometry (point) dofmap. For a cell, it /// gives the position in the point array of each local geometry node /// @param[in] element The element that describes the cell geometry map - /// @param[in] x The point coordinates. It is a `std::vector` - /// and uses row-major storage. The shape is `(num_points, 3)`. + /// @param[in] x The point coordinates. It is a `std::vector` and + /// uses row-major storage. The shape is `(num_points, 3)`. /// @param[in] dim The geometric dimension (`0 < dim <= 3`) /// @param[in] input_global_indices The 'global' input index of each /// point, commonly from a mesh input file. The type is @@ -71,6 +71,15 @@ class Geometry /// Move Assignment Geometry& operator=(Geometry&&) = default; + /// Copy constructor + template + Geometry astype() const + { + return Geometry(_index_map, _dofmap, _cmap, + std::vector(_x.begin(), _x.end()), _dim, + _input_global_indices); + } + /// Return Euclidean dimension of coordinate system int dim() const { return _dim; } @@ -131,8 +140,8 @@ class Geometry /// @brief Build Geometry from input data. /// /// This function should be called after the mesh topology is built. It -/// distributes the 'node' coordinate data to the required MPI process and then -/// creates a mesh::Geometry object. +/// distributes the 'node' coordinate data to the required MPI process +/// and then creates a mesh::Geometry object. /// /// @param[in] comm The MPI communicator to build the Geometry on /// @param[in] topology The mesh topology From b9576f62ea0048b822a8627712074e9ea340659b Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sat, 19 Nov 2022 19:27:21 +0000 Subject: [PATCH 03/11] Work on geometry type --- cpp/dolfinx/fem/assemble_matrix_impl.h | 56 +++++++++++++++++++------- 1 file changed, 42 insertions(+), 14 deletions(-) diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h index 9262d8068dd..72a848c9378 100644 --- a/cpp/dolfinx/fem/assemble_matrix_impl.h +++ b/cpp/dolfinx/fem/assemble_matrix_impl.h @@ -28,7 +28,7 @@ namespace dolfinx::fem::impl template void assemble_cells( la::MatSet auto mat_set, - const mesh::Geometry> geometry, + const mesh::Geometry> geometry, std::span cells, const std::function&, const std::span&, @@ -124,7 +124,7 @@ void assemble_cells( template void assemble_exterior_facets( la::MatSet auto mat_set, - const mesh::Geometry>& geometry, + const mesh::Geometry>& geometry, std::span facets, const std::function&, const std::span&, @@ -219,7 +219,7 @@ void assemble_exterior_facets( template void assemble_interior_facets( la::MatSet auto mat_set, - const mesh::Geometry>& geometry, int num_cell_facets, + const mesh::Geometry>& geometry, int num_cell_facets, std::span facets, const std::function&, const std::span&, @@ -363,7 +363,9 @@ void assemble_interior_facets( /// are applied. Matrix is not finalised. template void assemble_matrix( - la::MatSet auto mat_set, const Form& a, std::span constants, + la::MatSet auto mat_set, const Form& a, + const mesh::Geometry>& geometry, + std::span constants, const std::map, std::pair, int>>& coefficients, std::span bc0, std::span bc1) @@ -412,9 +414,9 @@ void assemble_matrix( const auto& fn = a.kernel(IntegralType::cell, i); const auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i}); const std::vector& cells = a.cell_domains(i); - impl::assemble_cells(mat_set, mesh->geometry(), cells, dof_transform, dofs0, - bs0, dof_transform_to_transpose, dofs1, bs1, bc0, bc1, - fn, coeffs, cstride, constants, cell_info); + impl::assemble_cells(mat_set, geometry, cells, dof_transform, dofs0, bs0, + dof_transform_to_transpose, dofs1, bs1, bc0, bc1, fn, + coeffs, cstride, constants, cell_info); } for (int i : a.integral_ids(IntegralType::exterior_facet)) @@ -423,10 +425,10 @@ void assemble_matrix( const auto& [coeffs, cstride] = coefficients.at({IntegralType::exterior_facet, i}); const std::vector& facets = a.exterior_facet_domains(i); - impl::assemble_exterior_facets( - mat_set, mesh->geometry(), facets, dof_transform, dofs0, bs0, - dof_transform_to_transpose, dofs1, bs1, bc0, bc1, fn, coeffs, cstride, - constants, cell_info); + impl::assemble_exterior_facets(mat_set, geometry, facets, dof_transform, + dofs0, bs0, dof_transform_to_transpose, + dofs1, bs1, bc0, bc1, fn, coeffs, cstride, + constants, cell_info); } if (a.num_integrals(IntegralType::interior_facet) > 0) @@ -452,11 +454,37 @@ void assemble_matrix( = coefficients.at({IntegralType::interior_facet, i}); const std::vector& facets = a.interior_facet_domains(i); impl::assemble_interior_facets( - mat_set, mesh->geometry(), num_cell_facets, facets, dof_transform, - *dofmap0, bs0, dof_transform_to_transpose, *dofmap1, bs1, bc0, bc1, - fn, coeffs, cstride, c_offsets, constants, cell_info, get_perm); + mat_set, geometry, num_cell_facets, facets, dof_transform, *dofmap0, + bs0, dof_transform_to_transpose, *dofmap1, bs1, bc0, bc1, fn, coeffs, + cstride, c_offsets, constants, cell_info, get_perm); } } } +/// The matrix A must already be initialised. The matrix may be a proxy, +/// i.e. a view into a larger matrix, and assembly is performed using +/// local indices. Rows (bc0) and columns (bc1) with Dirichlet +/// conditions are zeroed. Markers (bc0 and bc1) can be empty if not bcs +/// are applied. Matrix is not finalised. +template +void assemble_matrix( + la::MatSet auto mat_set, const Form& a, std::span constants, + const std::map, + std::pair, int>>& coefficients, + std::span bc0, std::span bc1) +{ + std::shared_ptr mesh = a.mesh(); + assert(mesh); + if constexpr (std::is_same_v>) + return assemble_matrix(mat_set, a, mesh->geometry(), constants, + coefficients, bc0, bc1); + else + { + + return assemble_matrix(mat_set, a, + mesh->geometry().astype>(), + constants, coefficients, bc0, bc1); + } +} + } // namespace dolfinx::fem::impl From 57369ca4ac3810974e0226b12a2c7bd59ba8c596 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sat, 19 Nov 2022 19:28:36 +0000 Subject: [PATCH 04/11] Formatting fixes --- cpp/dolfinx/fem/assemble_matrix_impl.h | 13 +++++++------ cpp/dolfinx/fem/assemble_scalar_impl.h | 1 - cpp/dolfinx/fem/assemble_vector_impl.h | 8 +++----- 3 files changed, 10 insertions(+), 12 deletions(-) diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h index 72a848c9378..7c47136c7a9 100644 --- a/cpp/dolfinx/fem/assemble_matrix_impl.h +++ b/cpp/dolfinx/fem/assemble_matrix_impl.h @@ -476,14 +476,15 @@ void assemble_matrix( std::shared_ptr mesh = a.mesh(); assert(mesh); if constexpr (std::is_same_v>) - return assemble_matrix(mat_set, a, mesh->geometry(), constants, - coefficients, bc0, bc1); + { + assemble_matrix(mat_set, a, mesh->geometry(), constants, coefficients, bc0, + bc1); + } else { - - return assemble_matrix(mat_set, a, - mesh->geometry().astype>(), - constants, coefficients, bc0, bc1); + assemble_matrix(mat_set, a, + mesh->geometry().astype>(), + constants, coefficients, bc0, bc1); } } diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index 21a2b3bbec8..b51a0356a5f 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -230,7 +230,6 @@ T assemble_scalar( return assemble_scalar(M, mesh->geometry(), constants, coefficients); else { - return assemble_scalar(M, mesh->geometry().astype>(), constants, coefficients); } diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index 83e79ee3e9a..eedceac82bb 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -1073,13 +1073,11 @@ void assemble_vector( std::shared_ptr mesh = L.mesh(); assert(mesh); if constexpr (std::is_same_v>) - return assemble_vector(b, L, mesh->geometry(), constants, coefficients); + assemble_vector(b, L, mesh->geometry(), constants, coefficients); else { - - return assemble_vector(b, L, - mesh->geometry().astype>(), - constants, coefficients); + assemble_vector(b, L, mesh->geometry().astype>(), + constants, coefficients); } } } // namespace dolfinx::fem::impl From cf63bd547c275b24d9aa2e34819b3adbddcc4ce8 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sat, 19 Nov 2022 19:58:02 +0000 Subject: [PATCH 05/11] Refactoring --- cpp/dolfinx/fem/assemble_matrix_impl.h | 27 -------- cpp/dolfinx/fem/assemble_scalar_impl.h | 18 ----- cpp/dolfinx/fem/assemble_vector_impl.h | 68 +++++++++---------- cpp/dolfinx/fem/assembler.h | 94 ++++++++++++++++++-------- 4 files changed, 95 insertions(+), 112 deletions(-) diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h index 7c47136c7a9..b4fe30ceb40 100644 --- a/cpp/dolfinx/fem/assemble_matrix_impl.h +++ b/cpp/dolfinx/fem/assemble_matrix_impl.h @@ -461,31 +461,4 @@ void assemble_matrix( } } -/// The matrix A must already be initialised. The matrix may be a proxy, -/// i.e. a view into a larger matrix, and assembly is performed using -/// local indices. Rows (bc0) and columns (bc1) with Dirichlet -/// conditions are zeroed. Markers (bc0 and bc1) can be empty if not bcs -/// are applied. Matrix is not finalised. -template -void assemble_matrix( - la::MatSet auto mat_set, const Form& a, std::span constants, - const std::map, - std::pair, int>>& coefficients, - std::span bc0, std::span bc1) -{ - std::shared_ptr mesh = a.mesh(); - assert(mesh); - if constexpr (std::is_same_v>) - { - assemble_matrix(mat_set, a, mesh->geometry(), constants, coefficients, bc0, - bc1); - } - else - { - assemble_matrix(mat_set, a, - mesh->geometry().astype>(), - constants, coefficients, bc0, bc1); - } -} - } // namespace dolfinx::fem::impl diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index b51a0356a5f..900056e7989 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -217,22 +217,4 @@ T assemble_scalar( return value; } -/// Assemble functional into an scalar -template -T assemble_scalar( - const fem::Form& M, std::span constants, - const std::map, - std::pair, int>>& coefficients) -{ - std::shared_ptr mesh = M.mesh(); - assert(mesh); - if constexpr (std::is_same_v>) - return assemble_scalar(M, mesh->geometry(), constants, coefficients); - else - { - return assemble_scalar(M, mesh->geometry().astype>(), - constants, coefficients); - } -} - } // namespace dolfinx::fem::impl diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index eedceac82bb..574d0d81dcd 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -37,7 +37,7 @@ namespace dolfinx::fem::impl /// @tparam _bs1 The block size of the trial function dof map. template void _lift_bc_cells( - std::span b, const mesh::Geometry>& geometry, + std::span b, const mesh::Geometry>& geometry, FEkernel auto kernel, std::span cells, const std::function&, const std::span&, @@ -49,8 +49,7 @@ void _lift_bc_cells( const graph::AdjacencyList& dofmap1, int bs1, std::span constants, std::span coeffs, int cstride, std::span cell_info, std::span bc_values1, - std::span bc_markers1, std::span x0, - double scale) + std::span bc_markers1, std::span x0, T scale) { assert(_bs0 < 0 or _bs0 == bs0); assert(_bs1 < 0 or _bs1 == bs1); @@ -66,8 +65,6 @@ void _lift_bc_cells( // Data structures used in bc application std::vector> coordinate_dofs(3 * num_dofs_g); std::vector Ae, be; - const scalar_value_type_t _scale - = static_cast>(scale); for (std::size_t index = 0; index < cells.size(); ++index) { std::int32_t c = cells[index]; @@ -148,7 +145,7 @@ void _lift_bc_cells( // const T _x0 = 0.0; // be -= Ae.col(bs1 * j + k) * scale * (bc - _x0); for (int m = 0; m < num_rows; ++m) - be[m] -= Ae[m * num_cols + _bs1 * j + k] * _scale * (bc - _x0); + be[m] -= Ae[m * num_cols + _bs1 * j + k] * scale * (bc - _x0); } } } @@ -164,7 +161,7 @@ void _lift_bc_cells( const T _x0 = x0.empty() ? 0.0 : x0[jj]; // be -= Ae.col(bs1 * j + k) * scale * (bc - _x0); for (int m = 0; m < num_rows; ++m) - be[m] -= Ae[m * num_cols + bs1 * j + k] * _scale * (bc - _x0); + be[m] -= Ae[m * num_cols + bs1 * j + k] * scale * (bc - _x0); } } } @@ -194,7 +191,7 @@ void _lift_bc_cells( /// @tparam _bs1 The block size of the trial function dof map. template void _lift_bc_exterior_facets( - std::span b, const mesh::Geometry>& geometry, + std::span b, const mesh::Geometry>& geometry, FEkernel auto kernel, std::span facets, const std::function&, const std::span&, @@ -206,8 +203,7 @@ void _lift_bc_exterior_facets( const graph::AdjacencyList& dofmap1, int bs1, std::span constants, std::span coeffs, int cstride, std::span cell_info, std::span bc_values1, - std::span bc_markers1, std::span x0, - double scale) + std::span bc_markers1, std::span x0, T scale) { if (facets.empty()) return; @@ -221,8 +217,6 @@ void _lift_bc_exterior_facets( std::vector> coordinate_dofs(3 * num_dofs_g); std::vector Ae, be; assert(facets.size() % 2 == 0); - const scalar_value_type_t _scale - = static_cast>(scale); for (std::size_t index = 0; index < facets.size(); index += 2) { std::int32_t cell = facets[index]; @@ -283,7 +277,7 @@ void _lift_bc_exterior_facets( const T _x0 = x0.empty() ? 0.0 : x0[jj]; // be -= Ae.col(bs1 * j + k) * scale * (bc - _x0); for (int m = 0; m < num_rows; ++m) - be[m] -= Ae[m * num_cols + bs1 * j + k] * _scale * (bc - _x0); + be[m] -= Ae[m * num_cols + bs1 * j + k] * scale * (bc - _x0); } } } @@ -302,7 +296,7 @@ void _lift_bc_exterior_facets( /// @tparam _bs1 The block size of the trial function dof map. template void _lift_bc_interior_facets( - std::span b, const mesh::Geometry>& geometry, + std::span b, const mesh::Geometry>& geometry, int num_cell_facets, FEkernel auto kernel, std::span facets, const std::function&, @@ -317,7 +311,7 @@ void _lift_bc_interior_facets( std::span cell_info, const std::function& get_perm, std::span bc_values1, std::span bc_markers1, - std::span x0, double scale) + std::span x0, T scale) { if (facets.empty()) return; @@ -338,8 +332,6 @@ void _lift_bc_interior_facets( std::vector dmapjoint0, dmapjoint1; assert(facets.size() % 4 == 0); - const scalar_value_type_t _scale - = static_cast>(scale); for (std::size_t index = 0; index < facets.size(); index += 4) { std::array cells = {facets[index], facets[index + 2]}; @@ -444,7 +436,7 @@ void _lift_bc_interior_facets( const T _x0 = x0.empty() ? 0.0 : x0[jj]; // be -= Ae.col(bs1 * j + k) * scale * (bc - _x0); for (int m = 0; m < num_rows; ++m) - be[m] -= Ae[m * num_cols + bs1 * j + k] * _scale * (bc - _x0); + be[m] -= Ae[m * num_cols + bs1 * j + k] * scale * (bc - _x0); } } } @@ -463,8 +455,8 @@ void _lift_bc_interior_facets( // be -= Ae.col(offset + bs1 * j + k) * scale * (bc - x0[jj]); for (int m = 0; m < num_rows; ++m) { - be[m] -= Ae[m * num_cols + offset + bs1 * j + k] * _scale - * (bc - _x0); + be[m] + -= Ae[m * num_cols + offset + bs1 * j + k] * scale * (bc - _x0); } } } @@ -732,12 +724,14 @@ void assemble_interior_facets( /// solution' in a Newton method /// @param[in] scale Scaling to apply template -void lift_bc(std::span b, const Form& a, std::span constants, +void lift_bc(std::span b, const Form& a, + const mesh::Geometry>& geometry, + std::span constants, const std::map, std::pair, int>>& coefficients, std::span bc_values1, std::span bc_markers1, std::span x0, - double scale) + T scale) { std::shared_ptr mesh = a.mesh(); assert(mesh); @@ -785,22 +779,22 @@ void lift_bc(std::span b, const Form& a, std::span constants, const std::vector& cells = a.cell_domains(i); if (bs0 == 1 and bs1 == 1) { - _lift_bc_cells(b, mesh->geometry(), kernel, cells, dof_transform, + _lift_bc_cells(b, geometry, kernel, cells, dof_transform, dofmap0, bs0, dof_transform_to_transpose, dofmap1, bs1, constants, coeffs, cstride, cell_info, bc_values1, bc_markers1, x0, scale); } else if (bs0 == 3 and bs1 == 3) { - _lift_bc_cells(b, mesh->geometry(), kernel, cells, dof_transform, + _lift_bc_cells(b, geometry, kernel, cells, dof_transform, dofmap0, bs0, dof_transform_to_transpose, dofmap1, bs1, constants, coeffs, cstride, cell_info, bc_values1, bc_markers1, x0, scale); } else { - _lift_bc_cells(b, mesh->geometry(), kernel, cells, dof_transform, dofmap0, - bs0, dof_transform_to_transpose, dofmap1, bs1, constants, + _lift_bc_cells(b, geometry, kernel, cells, dof_transform, dofmap0, bs0, + dof_transform_to_transpose, dofmap1, bs1, constants, coeffs, cstride, cell_info, bc_values1, bc_markers1, x0, scale); } @@ -812,7 +806,7 @@ void lift_bc(std::span b, const Form& a, std::span constants, const auto& [coeffs, cstride] = coefficients.at({IntegralType::exterior_facet, i}); const std::vector& facets = a.exterior_facet_domains(i); - _lift_bc_exterior_facets(b, mesh->geometry(), kernel, facets, dof_transform, + _lift_bc_exterior_facets(b, geometry, kernel, facets, dof_transform, dofmap0, bs0, dof_transform_to_transpose, dofmap1, bs1, constants, coeffs, cstride, cell_info, bc_values1, bc_markers1, x0, scale); @@ -840,11 +834,10 @@ void lift_bc(std::span b, const Form& a, std::span constants, const auto& [coeffs, cstride] = coefficients.at({IntegralType::interior_facet, i}); const std::vector& facets = a.interior_facet_domains(i); - _lift_bc_interior_facets(b, mesh->geometry(), num_cell_facets, kernel, - facets, dof_transform, dofmap0, bs0, - dof_transform_to_transpose, dofmap1, bs1, - constants, coeffs, cstride, cell_info, get_perm, - bc_values1, bc_markers1, x0, scale); + _lift_bc_interior_facets( + b, geometry, num_cell_facets, kernel, facets, dof_transform, dofmap0, + bs0, dof_transform_to_transpose, dofmap1, bs1, constants, coeffs, + cstride, cell_info, get_perm, bc_values1, bc_markers1, x0, scale); } } } @@ -871,11 +864,12 @@ void lift_bc(std::span b, const Form& a, std::span constants, template void apply_lifting( std::span b, const std::vector>> a, + const mesh::Geometry>& geometry, const std::vector>& constants, const std::vector, std::pair, int>>>& coeffs, const std::vector>>>& bcs1, - const std::vector>& x0, double scale) + const std::vector>& x0, T scale) { // FIXME: make changes to reactivate this check if (!x0.empty() and x0.size() != a.size()) @@ -914,13 +908,13 @@ void apply_lifting( if (!x0.empty()) { - lift_bc(b, *a[j], constants[j], coeffs[j], bc_values1, bc_markers1, - x0[j], scale); + lift_bc(b, *a[j], geometry, constants[j], coeffs[j], bc_values1, + bc_markers1, x0[j], scale); } else { - lift_bc(b, *a[j], constants[j], coeffs[j], bc_values1, bc_markers1, - std::span(), scale); + lift_bc(b, *a[j], geometry, constants[j], coeffs[j], bc_values1, + bc_markers1, std::span(), scale); } } } diff --git a/cpp/dolfinx/fem/assembler.h b/cpp/dolfinx/fem/assembler.h index 4500f5e0c8e..6edb0f64432 100644 --- a/cpp/dolfinx/fem/assembler.h +++ b/cpp/dolfinx/fem/assembler.h @@ -9,6 +9,7 @@ #include "assemble_matrix_impl.h" #include "assemble_scalar_impl.h" #include "assemble_vector_impl.h" +#include "utils.h" #include #include #include @@ -58,7 +59,16 @@ T assemble_scalar( const std::map, std::pair, int>>& coefficients) { - return impl::assemble_scalar(M, constants, coefficients); + std::shared_ptr mesh = M.mesh(); + assert(mesh); + if constexpr (std::is_same_v>) + return impl::assemble_scalar(M, mesh->geometry(), constants, coefficients); + else + { + return impl::assemble_scalar( + M, mesh->geometry().astype>(), constants, + coefficients); + } } /// Assemble functional into scalar @@ -136,9 +146,21 @@ void apply_lifting( const std::vector, std::pair, int>>>& coeffs, const std::vector>>>& bcs1, - const std::vector>& x0, double scale) + const std::vector>& x0, T scale) { - impl::apply_lifting(b, a, constants, coeffs, bcs1, x0, scale); + std::shared_ptr mesh = a[0]->mesh(); + assert(mesh); + if constexpr (std::is_same_v>) + { + impl::apply_lifting(b, a, mesh->geometry(), constants, coeffs, bcs1, x0, + scale); + } + else + { + impl::apply_lifting( + b, a, mesh->geometry().astype>(), + constants, coeffs, bcs1, x0, scale); + } } /// Modify b such that: @@ -186,11 +208,48 @@ void apply_lifting( _coeffs; std::transform(coeffs.cbegin(), coeffs.cend(), std::back_inserter(_coeffs), [](auto& c) { return make_coefficients_span(c); }); + apply_lifting(b, a, _constants, _coeffs, bcs1, x0, scale); } // -- Matrices --------------------------------------------------------------- +/// @brief Assemble bilinear form into a matrix. Matrix must already be +/// initialised. Does not zero or finalise the matrix. +/// @param[in] mat_add The function for adding values into the matrix +/// @param[in] a The bilinear form to assemble +/// @param[in] constants Constants that appear in `a` +/// @param[in] coefficients Coefficients that appear in `a` +/// @param[in] dof_marker0 Boundary condition markers for the rows. If +/// bc[i] is true then rows i in A will be zeroed. The index i is a +/// local index. +/// @param[in] dof_marker1 Boundary condition markers for the columns. +/// If bc[i] is true then rows i in A will be zeroed. The index i is a +/// local index. +template +void assemble_matrix( + auto mat_add, const Form& a, std::span constants, + const std::map, + std::pair, int>>& coefficients, + std::span dof_marker0, + std::span dof_marker1) + +{ + std::shared_ptr mesh = a.mesh(); + assert(mesh); + if constexpr (std::is_same_v>) + { + impl::assemble_matrix(mat_add, a, mesh->geometry(), constants, coefficients, + dof_marker0, dof_marker1); + } + else + { + impl::assemble_matrix( + mat_add, a, mesh->geometry().astype>(), + constants, coefficients, dof_marker0, dof_marker1); + } +} + /// Assemble bilinear form into a matrix /// @param[in] mat_add The function for adding values into the matrix /// @param[in] a The bilinear from to assemble @@ -235,8 +294,8 @@ void assemble_matrix( } // Assemble - impl::assemble_matrix(mat_add, a, constants, coefficients, dof_marker0, - dof_marker1); + assemble_matrix(mat_add, a, constants, coefficients, dof_marker0, + dof_marker1); } /// Assemble bilinear form into a matrix @@ -259,31 +318,6 @@ void assemble_matrix( make_coefficients_span(coefficients), bcs); } -/// @brief Assemble bilinear form into a matrix. Matrix must already be -/// initialised. Does not zero or finalise the matrix. -/// @param[in] mat_add The function for adding values into the matrix -/// @param[in] a The bilinear form to assemble -/// @param[in] constants Constants that appear in `a` -/// @param[in] coefficients Coefficients that appear in `a` -/// @param[in] dof_marker0 Boundary condition markers for the rows. If -/// bc[i] is true then rows i in A will be zeroed. The index i is a -/// local index. -/// @param[in] dof_marker1 Boundary condition markers for the columns. -/// If bc[i] is true then rows i in A will be zeroed. The index i is a -/// local index. -template -void assemble_matrix( - auto mat_add, const Form& a, std::span constants, - const std::map, - std::pair, int>>& coefficients, - std::span dof_marker0, - std::span dof_marker1) - -{ - impl::assemble_matrix(mat_add, a, constants, coefficients, dof_marker0, - dof_marker1); -} - /// @brief Assemble bilinear form into a matrix. Matrix must already be /// initialised. Does not zero or finalise the matrix. /// @param[in] mat_add The function for adding values into the matrix From 2be03d6d1241037b41a341284ff42fecf90206d1 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sat, 19 Nov 2022 21:10:56 +0000 Subject: [PATCH 06/11] Doc fixes --- cpp/dolfinx/fem/assemble_vector_impl.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index 574d0d81dcd..2af569f7e68 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -715,6 +715,7 @@ void assemble_interior_facets( /// /// @param[in,out] b The vector to be modified /// @param[in] a The bilinear form that generates A +/// @param[in] geometry The mesh geometry /// @param[in] constants Constants that appear in `a` /// @param[in] coefficients Coefficients that appear in `a` /// @param[in] bc_values1 The boundary condition 'values' @@ -854,6 +855,7 @@ void lift_bc(std::span b, const Form& a, /// @param[in,out] b The vector to be modified /// @param[in] a The bilinear forms, where a[j] is the form that /// generates A_j +/// @param[in] geometry The mesh geometry /// @param[in] constants Constants that appear in `a` /// @param[in] coeffs Coefficients that appear in `a` /// @param[in] bcs1 List of boundary conditions for each block, i.e. @@ -924,6 +926,7 @@ void apply_lifting( /// @param[in,out] b The vector to be assembled. It will not be zeroed /// before assembly. /// @param[in] L The linear forms to assemble into b +/// @param[in] geometry The mesh geometry /// @param[in] constants Packed constants that appear in `L` /// @param[in] coefficients Packed coefficients that appear in `L` template From 6f43bb04a6bf05430a9d7bb00de931e9d31ef2fa Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sun, 20 Nov 2022 16:14:46 +0000 Subject: [PATCH 07/11] Minor inteface update --- cpp/dolfinx/fem/assembler.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cpp/dolfinx/fem/assembler.h b/cpp/dolfinx/fem/assembler.h index 6edb0f64432..440bdaa3388 100644 --- a/cpp/dolfinx/fem/assembler.h +++ b/cpp/dolfinx/fem/assembler.h @@ -179,7 +179,7 @@ template void apply_lifting( std::span b, const std::vector>>& a, const std::vector>>>& bcs1, - const std::vector>& x0, double scale) + const std::vector>& x0, T scale) { std::vector< std::map, std::pair, int>>> @@ -412,7 +412,7 @@ void set_diagonal(auto set_fn, const fem::FunctionSpace& V, template void set_bc(std::span b, const std::vector>>& bcs, - std::span x0, double scale = 1.0) + std::span x0, T scale = 1) { if (b.size() > x0.size()) throw std::runtime_error("Size mismatch between b and x0 vectors."); @@ -429,7 +429,7 @@ void set_bc(std::span b, template void set_bc(std::span b, const std::vector>>& bcs, - double scale = 1.0) + T scale = 1) { for (const auto& bc : bcs) { From 6a4f9a97173ad42824abd49400771ee6408f537c Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sun, 20 Nov 2022 16:23:20 +0000 Subject: [PATCH 08/11] Small type updates --- cpp/dolfinx/fem/DirichletBC.h | 2 +- cpp/dolfinx/fem/petsc.cpp | 6 +++--- cpp/dolfinx/fem/petsc.h | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/cpp/dolfinx/fem/DirichletBC.h b/cpp/dolfinx/fem/DirichletBC.h index 9e16a4f9e29..06437a570cd 100644 --- a/cpp/dolfinx/fem/DirichletBC.h +++ b/cpp/dolfinx/fem/DirichletBC.h @@ -348,7 +348,7 @@ class DirichletBC /// of the array @p x should be equal to the number of dofs owned by /// this rank. /// @param[in] scale The scaling value to apply - void set(std::span x, double scale = 1.0) const + void set(std::span x, T scale = 1) const { if (std::holds_alternative>>(_g)) { diff --git a/cpp/dolfinx/fem/petsc.cpp b/cpp/dolfinx/fem/petsc.cpp index 6cf4d0a986a..74426c04820 100644 --- a/cpp/dolfinx/fem/petsc.cpp +++ b/cpp/dolfinx/fem/petsc.cpp @@ -304,7 +304,7 @@ void fem::petsc::apply_lifting( coeffs, const std::vector< std::vector>>>& bcs1, - const std::vector& x0, double scale) + const std::vector& x0, PetscScalar scale) { Vec b_local; VecGhostGetLocalForm(b, &b_local); @@ -350,7 +350,7 @@ void fem::petsc::apply_lifting( Vec b, const std::vector>>& a, const std::vector< std::vector>>>& bcs1, - const std::vector& x0, double scale) + const std::vector& x0, PetscScalar scale) { Vec b_local; VecGhostGetLocalForm(b, &b_local); @@ -394,7 +394,7 @@ void fem::petsc::apply_lifting( void fem::petsc::set_bc( Vec b, const std::vector>>& bcs, - const Vec x0, double scale) + const Vec x0, PetscScalar scale) { PetscInt n = 0; VecGetLocalSize(b, &n); diff --git a/cpp/dolfinx/fem/petsc.h b/cpp/dolfinx/fem/petsc.h index f6bc374592f..a966b5c722d 100644 --- a/cpp/dolfinx/fem/petsc.h +++ b/cpp/dolfinx/fem/petsc.h @@ -163,6 +163,6 @@ void apply_lifting( void set_bc( Vec b, const std::vector>>& bcs, - const Vec x0, double scale = 1.0); + const Vec x0, PetscScalar scale = 1); } // namespace petsc } // namespace dolfinx::fem From 98d357a84c2d8754440c698fa39ad54ff5db7b34 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sun, 20 Nov 2022 18:01:31 +0000 Subject: [PATCH 09/11] Type fixes --- cpp/dolfinx/fem/DirichletBC.h | 2 +- cpp/dolfinx/fem/petsc.h | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/cpp/dolfinx/fem/DirichletBC.h b/cpp/dolfinx/fem/DirichletBC.h index 06437a570cd..88f9fc7f200 100644 --- a/cpp/dolfinx/fem/DirichletBC.h +++ b/cpp/dolfinx/fem/DirichletBC.h @@ -385,7 +385,7 @@ class DirichletBC /// @param[in] x The array in which to set `scale * (x0 - x_bc)` /// @param[in] x0 The array used in compute the value to set /// @param[in] scale The scaling value to apply - void set(std::span x, std::span x0, double scale = 1.0) const + void set(std::span x, std::span x0, T scale = 1) const { if (std::holds_alternative>>(_g)) { diff --git a/cpp/dolfinx/fem/petsc.h b/cpp/dolfinx/fem/petsc.h index a966b5c722d..ff65c243cc5 100644 --- a/cpp/dolfinx/fem/petsc.h +++ b/cpp/dolfinx/fem/petsc.h @@ -124,7 +124,7 @@ void apply_lifting( coeffs, const std::vector< std::vector>>>& bcs1, - const std::vector& x0, double scale); + const std::vector& x0, PetscScalar scale); // FIXME: clarify how x0 is used // FIXME: if bcs entries are set @@ -148,7 +148,7 @@ void apply_lifting( Vec b, const std::vector>>& a, const std::vector< std::vector>>>& bcs1, - const std::vector& x0, double scale); + const std::vector& x0, PetscScalar scale); // -- Setting bcs ------------------------------------------------------------ From fe8b21ed77a0d7ec740f30d164723885d78e171f Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sun, 20 Nov 2022 18:15:27 +0000 Subject: [PATCH 10/11] Demo updates --- cpp/demo/poisson/main.cpp | 2 +- cpp/demo/poisson_matrix_free/main.cpp | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/cpp/demo/poisson/main.cpp b/cpp/demo/poisson/main.cpp index a4532d1cd43..efe72dd1dae 100644 --- a/cpp/demo/poisson/main.cpp +++ b/cpp/demo/poisson/main.cpp @@ -229,7 +229,7 @@ int main(int argc, char* argv[]) b.set(0.0); fem::assemble_vector(b.mutable_array(), *L); - fem::apply_lifting(b.mutable_array(), {a}, {{bc}}, {}, 1.0); + fem::apply_lifting(b.mutable_array(), {a}, {{bc}}, {}, T(1)); b.scatter_rev(std::plus()); fem::set_bc(b.mutable_array(), {bc}); diff --git a/cpp/demo/poisson_matrix_free/main.cpp b/cpp/demo/poisson_matrix_free/main.cpp index 004dfee2732..0e1f2a05707 100644 --- a/cpp/demo/poisson_matrix_free/main.cpp +++ b/cpp/demo/poisson_matrix_free/main.cpp @@ -164,14 +164,14 @@ int main(int argc, char* argv[]) // Apply lifting to account for Dirichlet boundary condition // b <- b - A * x_bc - fem::set_bc(ui->x()->mutable_array(), {bc}, -1.0); + fem::set_bc(ui->x()->mutable_array(), {bc}, T(-1)); fem::assemble_vector(b.mutable_array(), *M); // Communicate ghost values b.scatter_rev(std::plus()); // Set BC dofs to zero (effectively zeroes columns of A) - fem::set_bc(b.mutable_array(), {bc}, 0.0); + fem::set_bc(b.mutable_array(), {bc}, T(0)); b.scatter_fwd(); @@ -196,7 +196,7 @@ int main(int argc, char* argv[]) fem::make_coefficients_span(coeff)); // Set BC dofs to zero (effectively zeroes rows of A) - fem::set_bc(y.mutable_array(), {bc}, 0.0); + fem::set_bc(y.mutable_array(), {bc}, T(0)); // Accumuate ghost values y.scatter_rev(std::plus()); @@ -210,7 +210,7 @@ int main(int argc, char* argv[]) int num_it = linalg::cg(*u->x(), b, action, 200, 1e-6); // Set BC values in the solution vectors - fem::set_bc(u->x()->mutable_array(), {bc}, 1.0); + fem::set_bc(u->x()->mutable_array(), {bc}, T(1)); // Compute L2 error (squared) of the solution vector e = (u - u_d, u // - u_d)*dx From 46ef06cc2ef2c3d5c039c15d376cfbfd74c9b48d Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Wed, 1 Feb 2023 07:57:07 +0000 Subject: [PATCH 11/11] Remove duplication --- cpp/dolfinx/fem/assemble_scalar_impl.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index 900056e7989..a72db267ead 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -24,7 +24,7 @@ namespace dolfinx::fem::impl /// Assemble functional over cells template T assemble_cells( - const mesh::Geometry>>& geometry, + const mesh::Geometry>& geometry, std::span cells, FEkernel auto fn, std::span constants, std::span coeffs, int cstride) {