diff --git a/cpp/demo/hyperelasticity/main.cpp b/cpp/demo/hyperelasticity/main.cpp index a9458bfbc1a..4e55588dade 100644 --- a/cpp/demo/hyperelasticity/main.cpp +++ b/cpp/demo/hyperelasticity/main.cpp @@ -20,10 +20,10 @@ class HyperElasticProblem : public nls::NonlinearProblem std::shared_ptr> J, std::vector>> bcs) : _u(u), _l(L), _j(J), _bcs(bcs), - _b(L->function_space(0)->dofmap()->index_map), + _b(L->function_spaces()[0]->dofmap()->index_map), _matA(fem::create_matrix(*J)) { - auto map = L->function_space(0)->dofmap()->index_map; + auto map = L->function_spaces()[0]->dofmap()->index_map; const int bs = map->block_size(); std::int32_t size_local = bs * map->size_local(); std::int32_t num_ghosts = bs * map->num_ghosts(); @@ -82,7 +82,7 @@ class HyperElasticProblem : public nls::NonlinearProblem MatZeroEntries(_matA.mat()); fem::assemble_matrix(la::PETScMatrix::add_fn(_matA.mat()), *_j, _bcs); fem::add_diagonal(la::PETScMatrix::add_fn(_matA.mat()), - *_j->function_space(0), _bcs); + *_j->function_spaces()[0], _bcs); _matA.apply(la::PETScMatrix::AssemblyType::FINAL); return _matA.mat(); } @@ -125,9 +125,10 @@ int main(int argc, char* argv[]) // Define solution function auto u = std::make_shared>(V); - auto a - = fem::create_form(create_form_hyperelasticity_J, {V, V}); - auto L = fem::create_form(create_form_hyperelasticity_F, {V}); + auto a = fem::create_form(create_form_hyperelasticity_J, + {V, V}, {{"u", u}}, {}, {}); + auto L = fem::create_form(create_form_hyperelasticity_F, {V}, + {{"u", u}}, {}, {}); auto u_rotation = std::make_shared>(V); u_rotation->interpolate([](auto& x) { @@ -165,9 +166,6 @@ int main(int argc, char* argv[]) Eigen::RowMajor>::Zero(3, x.cols()); }); - L->set_coefficients({{"u", u}}); - a->set_coefficients({{"u", u}}); - // Create Dirichlet boundary conditions auto u0 = std::make_shared>(V); diff --git a/cpp/demo/poisson/main.cpp b/cpp/demo/poisson/main.cpp index 638a447252e..d19c1e8327d 100644 --- a/cpp/demo/poisson/main.cpp +++ b/cpp/demo/poisson/main.cpp @@ -134,13 +134,18 @@ int main(int argc, char* argv[]) // // .. code-block:: cpp - // Define variational forms - auto a = fem::create_form(create_form_poisson_a, {V, V}); - auto L = fem::create_form(create_form_poisson_L, {V}); + // Prepare and set Constants for the bilinear form + auto kappa = std::make_shared>(2.0); auto f = std::make_shared>(V); auto g = std::make_shared>(V); + // Define variational forms + auto a = fem::create_form(create_form_poisson_a, {V, V}, {}, + {{"kappa", kappa}}, {}); + auto L = fem::create_form(create_form_poisson_L, {V}, + {{"f", f}, {"g", g}}, {}, {}); + // Now, the Dirichlet boundary condition (:math:`u = 0`) can be created // using the class :cpp:class:`DirichletBC`. A :cpp:class:`DirichletBC` // takes two arguments: the value of the boundary condition, @@ -173,11 +178,6 @@ int main(int argc, char* argv[]) }); g->interpolate([](auto& x) { return Eigen::sin(5 * x.row(0)); }); - L->set_coefficients({{"f", f}, {"g", g}}); - - // Prepare and set Constants for the bilinear form - auto kappa = std::make_shared>(2.0); - a->set_constants({{"kappa", kappa}}); // Now, we have specified the variational forms and can consider the // solution of the variational problem. First, we need to define a @@ -191,7 +191,7 @@ int main(int argc, char* argv[]) // Compute solution function::Function u(V); la::PETScMatrix A = fem::create_matrix(*a); - la::PETScVector b(*L->function_space(0)->dofmap()->index_map); + la::PETScVector b(*L->function_spaces()[0]->dofmap()->index_map); MatZeroEntries(A.mat()); fem::assemble_matrix(la::PETScMatrix::add_fn(A.mat()), *a, bc); diff --git a/cpp/dolfinx/fem/CMakeLists.txt b/cpp/dolfinx/fem/CMakeLists.txt index f663afe7971..0b2d6ae3efe 100644 --- a/cpp/dolfinx/fem/CMakeLists.txt +++ b/cpp/dolfinx/fem/CMakeLists.txt @@ -15,8 +15,6 @@ set(HEADERS_fem ${CMAKE_CURRENT_SOURCE_DIR}/utils.h ${CMAKE_CURRENT_SOURCE_DIR}/FiniteElement.h ${CMAKE_CURRENT_SOURCE_DIR}/Form.h - ${CMAKE_CURRENT_SOURCE_DIR}/FormCoefficients.h - ${CMAKE_CURRENT_SOURCE_DIR}/FormIntegrals.h ${CMAKE_CURRENT_SOURCE_DIR}/ReferenceCellGeometry.h ${CMAKE_CURRENT_SOURCE_DIR}/SparsityPatternBuilder.h PARENT_SCOPE) diff --git a/cpp/dolfinx/fem/FiniteElement.h b/cpp/dolfinx/fem/FiniteElement.h index 01362559595..cd6d21ba7ab 100644 --- a/cpp/dolfinx/fem/FiniteElement.h +++ b/cpp/dolfinx/fem/FiniteElement.h @@ -28,9 +28,21 @@ class FiniteElement /// @param[in] ufc_element UFC finite element explicit FiniteElement(const ufc_finite_element& ufc_element); + /// Copy constructor + FiniteElement(const FiniteElement& element) = default; + + /// Move constructor + FiniteElement(FiniteElement&& element) = default; + /// Destructor virtual ~FiniteElement() = default; + /// Copy assignment + FiniteElement& operator=(const FiniteElement& element) = default; + + /// Move assignment + FiniteElement& operator=(FiniteElement&& element) = default; + /// String identifying the finite element /// @return Element signature std::string signature() const; diff --git a/cpp/dolfinx/fem/Form.h b/cpp/dolfinx/fem/Form.h index d553bda2cdd..55d2f8da217 100644 --- a/cpp/dolfinx/fem/Form.h +++ b/cpp/dolfinx/fem/Form.h @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2014 Anders Logg +// Copyright (C) 2019-2020 Garth N. Wells and Chris Richardson // // This file is part of DOLFINX (https://www.fenicsproject.org) // @@ -6,14 +6,11 @@ #pragma once -#include "FormCoefficients.h" -#include "FormIntegrals.h" -#include -#include +#include +#include +#include #include -#include #include -#include #include #include @@ -26,18 +23,19 @@ template class Constant; template class Function; -class FunctionSpace; } // namespace function -namespace mesh +namespace fem { -class Mesh; -template -class MeshTags; -} // namespace mesh -namespace fem +/// Type of integral +enum class IntegralType : std::int8_t { + cell = 0, + exterior_facet = 1, + interior_facet = 2, + vertex = 3 +}; /// Class for variational forms /// @@ -69,49 +67,70 @@ class Form public: /// Create form /// - /// @param[in] function_spaces Function Spaces - /// @param[in] integrals + /// @param[in] function_spaces Function spaces for the form arguments + /// @param[in] integrals The integrals in the form. The first key is + /// the domain type. For each key there is a pair (list[domain id, + /// integration kernel], domain markers). /// @param[in] coefficients - /// @param[in] constants Vector of pairs (name, constant). The index - /// in the vector is the position of the constant in the original - /// (nonsimplified) form. + /// @param[in] constants Constants in the Form + /// @param[in] needs_permutation_data Set to true is any of the + /// integration kernels require cell permutation data + /// @param[in] mesh The mesh of the domain. This is required when + /// there are not argument functions from which the mesh can be + /// extracted, e.g. for functionals Form(const std::vector>& function_spaces, - const FormIntegrals& integrals, - const FormCoefficients& coefficients, - const std::vector< - std::pair>>> - constants) - : _integrals(integrals), _coefficients(coefficients), - _constants(constants), _function_spaces(function_spaces) + const std::map< + IntegralType, + std::pair< + std::vector>>, + const mesh::MeshTags*>>& integrals, + const std::vector>>& + coefficients, + const std::vector>>& + constants, + bool needs_permutation_data, + const std::shared_ptr& mesh = nullptr) + : _function_spaces(function_spaces), _coefficients(coefficients), + _constants(constants), _mesh(mesh), + _needs_permutation_data(needs_permutation_data) { - // Set _mesh from function::FunctionSpace, and check they are the same - if (!function_spaces.empty()) + // Extract _mesh from function::FunctionSpace, and check they are the same + if (!_mesh and !function_spaces.empty()) _mesh = function_spaces[0]->mesh(); for (const auto& V : function_spaces) if (_mesh != V->mesh()) throw std::runtime_error("Incompatible mesh"); + if (!_mesh) + throw std::runtime_error("No mesh could be associated with the Form."); - // Set markers for default integrals - if (_mesh) - _integrals.set_default_domains(*_mesh); - } + // Store kernels, looping over integrals by domain type (dimension) + for (auto& integral_type : integrals) + { + // Add key to map + const IntegralType type = integral_type.first; + auto it = _integrals.emplace( + type, std::map>>()); + + // Loop over integrals kernels + for (auto& integral : integral_type.second.first) + it.first->second.insert({integral.first, {integral.second, {}}}); + + // FIXME: do this neatly via a static function + // Set domains for integral type + if (integral_type.second.second) + { + assert(_mesh == integral_type.second.second->mesh()); + set_domains(type, *integral_type.second.second); + } + } - /// @warning Experimental - /// - /// Create form (no UFC integrals). Integrals can be attached later - /// using FormIntegrals::set_cell_tabulate_tensor. - /// - /// @param[in] function_spaces Vector of function spaces - /// @param[in] need_mesh_permutation_data Set to true if mesh entity - /// permutation data is required - Form(const std::vector>& - function_spaces, - bool need_mesh_permutation_data) - : Form(function_spaces, FormIntegrals({}, need_mesh_permutation_data), - FormCoefficients({}), {}) - { - // Do nothing + // FIXME: do this neatly via a static function + // Set markers for default integrals + set_default_domains(*_mesh); } /// Copy constructor @@ -128,153 +147,331 @@ class Form /// @return The rank of the form int rank() const { return _function_spaces.size(); } - /// Set coefficient with given name - /// @param[in] coefficients Map from coefficient name to the - /// coefficient - void set_coefficients( - const std::map>>& - coefficients) - { - std::for_each(coefficients.begin(), coefficients.end(), - [this](auto& c) { _coefficients.set(c.first, c.second); }); - } + /// Extract common mesh for the form + /// @return The mesh + std::shared_ptr mesh() const { return _mesh; } - /// Set constants based on their names. Names of the constants must - /// agree with their names in UFL file. - void set_constants( - const std::map>>& - constants) + /// Return function spaces for all arguments + /// @return Function spaces + const std::vector>& + function_spaces() const { - for (auto const& constant : constants) - { - // Find matching string in existing constants - const std::string& name = constant.first; - auto it - = std::find_if(_constants.begin(), _constants.end(), - [&name](const auto& q) { return q.first == name; }); - if (it != _constants.end()) - it->second = constant.second; - else - throw std::runtime_error("Constant '" + name + "' not found in form"); - } + return _function_spaces; } - /// Check if all constants associated with the form have been set - /// @return True if all Form constants have been set - bool all_constants_set() const + /// Get the function for 'kernel' for integral i of given + /// type + /// @param[in] type Integral type + /// @param[in] i Integral number + /// @return Function to call for tabulate_tensor + const std::function& + kernel(IntegralType type, int i) const { - for (const auto& constant : _constants) - if (!constant.second) - return false; - return true; + auto it0 = _integrals.find(type); + if (it0 == _integrals.end()) + throw std::runtime_error("No kernels for requested type."); + auto it1 = it0->second.find(i); + if (it1 == it0->second.end()) + throw std::runtime_error("No kernel for requested index."); + + return it1->second.first; } - /// Return names of any constants that have not been set - /// @return Names of unset constants - std::set get_unset_constants() const + /// Get types of integrals in the form + /// @return Integrals types + std::set integral_types() const { - std::set unset; - std::for_each(_constants.begin(), _constants.end(), [&unset](auto& c) { - if (!c.second) - unset.insert(c.first); - }); - return unset; + std::set set; + for (auto& type : _integrals) + set.insert(type.first); + return set; } - /// @todo Remove this function and make sure the mesh can be set via - /// the constructor - /// - /// Set mesh, necessary for functionals when there are no function - /// spaces - /// @param[in] mesh The mesh - void set_mesh(const std::shared_ptr& mesh) + /// Number of integrals of given type + /// @param[in] type Integral type + /// @return Number of integrals + int num_integrals(IntegralType type) const { - _mesh = mesh; - // Set markers for default integrals - _integrals.set_default_domains(*_mesh); + if (auto it = _integrals.find(type); it == _integrals.end()) + return 0; + else + return it->second.size(); } - /// Extract common mesh from form - /// @return The mesh - std::shared_ptr mesh() const { return _mesh; } - - /// Return function space for given argument - /// @param[in] i Index of the argument - /// @return Function space - std::shared_ptr function_space(int i) const + /// Get the IDs for integrals (kernels) for given integral type. The + /// IDs correspond to the domain IDs which the integrals are defined + /// for in the form. ID=-1 is the default integral over the whole + /// domain. + /// @param[in] type Integral type + /// @return List of IDs for given integral type + std::vector integral_ids(IntegralType type) const { - return _function_spaces.at(i); + std::vector ids; + if (auto it = _integrals.find(type); it != _integrals.end()) + { + for (auto& kernel : it->second) + ids.push_back(kernel.first); + } + return ids; } - /// Return function spaces for all arguments - /// @return Function spaces - std::vector> - function_spaces() const + /// Get the list of mesh entity indices for the ith integral (kernel) + /// for the given domain type, i.e. for cell integrals a list of cell + /// indices, for facet integrals a list of facet indices, etc. + /// @param[in] type The integral type + /// @param[in] i Integral (kernel) index + /// @return List of active entities for the given integral (kernel) + const std::vector& domains(IntegralType type, int i) const { - return _function_spaces; + auto it0 = _integrals.find(type); + if (it0 == _integrals.end()) + throw std::runtime_error("No kernels for requested type."); + auto it1 = it0->second.find(i); + if (it1 == it0->second.end()) + throw std::runtime_error("No kernel for requested index."); + return it1->second.second; } - /// Register the function for 'tabulate_tensor' for cell integral i - void set_tabulate_tensor( - IntegralType type, int i, - const std::function& fn) + /// Access coefficients + const std::vector>> + coefficients() const { - _integrals.set_tabulate_tensor(type, i, fn); - if (i == -1 and _mesh) - _integrals.set_default_domains(*_mesh); + return _coefficients; } - /// Access coefficients - FormCoefficients& coefficients() { return _coefficients; } - - /// Access coefficients - const FormCoefficients& coefficients() const { return _coefficients; } - - /// Access form integrals - const FormIntegrals& integrals() const { return _integrals; } + /// Get bool indicating whether permutation data needs to be passed + /// into these integrals + /// @return True if cell permutation data is required + bool needs_permutation_data() const { return _needs_permutation_data; } - /// Access constants - /// @return Vector of attached constants with their names. Names are - /// used to set constants in user's c++ code. Index in the vector is - /// the position of the constant in the original (nonsimplified) form. - std::vector< - std::pair>>>& - constants() + /// Offset for each coefficient expansion array on a cell. Used to + /// pack data for multiple coefficients in a flat array. The last + /// entry is the size required to store all coefficients. + std::vector coefficient_offsets() const { - return _constants; + std::vector n{0}; + for (const auto& c : _coefficients) + { + if (!c) + throw std::runtime_error("Not all form coefficients have been set."); + n.push_back(n.back() + c->function_space()->element()->space_dimension()); + } + return n; } /// Access constants - /// @return Vector of attached constants with their names. Names are - /// used to set constants in user's c++ code. Index in the vector is - /// the position of the constant in the original (nonsimplified) form. - const std::vector< - std::pair>>>& + const std::vector>>& constants() const { return _constants; } private: - // Integrals associated with the Form - FormIntegrals _integrals; + /// Sets the entity indices to assemble over for kernels with a domain ID. + /// @param[in] type Integral type + /// @param[in] marker MeshTags with domain ID. Entities with marker + /// 'i' will be assembled over using the kernel with ID 'i'. The + /// MeshTags is not stored. + void set_domains(IntegralType type, const mesh::MeshTags& marker) + { + auto it0 = _integrals.find(type); + assert(it0 != _integrals.end()); + + std::shared_ptr mesh = marker.mesh(); + const mesh::Topology& topology = mesh->topology(); + const int tdim = topology.dim(); + int dim = tdim; + if (type == IntegralType::exterior_facet + or type == IntegralType::interior_facet) + { + dim = tdim - 1; + mesh->topology_mutable().create_connectivity(dim, tdim); + } + else if (type == IntegralType::vertex) + dim = 0; - // Coefficients associated with the Form - FormCoefficients _coefficients; + if (dim != marker.dim()) + { + throw std::runtime_error("Invalid MeshTags dimension:" + + std::to_string(marker.dim())); + } - // Constants associated with the Form - std::vector< - std::pair>>> - _constants; + // Get all integrals for considered entity type + std::map>>& integrals + = it0->second; + + // Get mesh tag data + const std::vector& values = marker.values(); + const std::vector& tagged_entities = marker.indices(); + assert(topology.index_map(dim)); + const auto entity_end + = std::lower_bound(tagged_entities.begin(), tagged_entities.end(), + topology.index_map(dim)->size_local()); + + if (dim == tdim - 1) + { + auto f_to_c = topology.connectivity(tdim - 1, tdim); + assert(f_to_c); + if (type == IntegralType::exterior_facet) + { + // Only need to consider shared facets when there are no ghost + // cells + assert(topology.index_map(tdim)); + std::set fwd_shared; + if (topology.index_map(tdim)->num_ghosts() == 0) + { + fwd_shared.insert( + topology.index_map(tdim - 1)->shared_indices().begin(), + topology.index_map(tdim - 1)->shared_indices().end()); + } + + for (auto f = tagged_entities.begin(); f != entity_end; ++f) + { + // All "owned" facets connected to one cell, that are not + // shared, should be external + if (f_to_c->num_links(*f) == 1 + and fwd_shared.find(*f) == fwd_shared.end()) + { + const std::size_t i = std::distance(tagged_entities.cbegin(), f); + if (auto it = integrals.find(values[i]); it != integrals.end()) + it->second.second.push_back(*f); + } + } + } + else if (type == IntegralType::interior_facet) + { + for (auto f = tagged_entities.begin(); f != entity_end; ++f) + { + if (f_to_c->num_links(*f) == 2) + { + const std::size_t i = std::distance(tagged_entities.cbegin(), f); + if (auto it = integrals.find(values[i]); it != integrals.end()) + it->second.second.push_back(*f); + } + } + } + } + else + { + // For cell and vertex integrals use all markers (but not on ghost + // entities) + for (auto e = tagged_entities.begin(); e != entity_end; ++e) + { + const std::size_t i = std::distance(tagged_entities.cbegin(), e); + if (auto it = integrals.find(values[i]); it != integrals.end()) + it->second.second.push_back(*e); + } + } + } + + /// If there exists a default integral of any type, set the list of + /// entities for those integrals from the mesh topology. For cell + /// integrals, this is all cells. For facet integrals, it is either + /// all interior or all exterior facets. + /// @param[in] mesh Mesh + void set_default_domains(const mesh::Mesh& mesh) + { + const mesh::Topology& topology = mesh.topology(); + const int tdim = topology.dim(); + + // Cells. If there is a default integral, define it on all owned cells + if (auto kernels = _integrals.find(IntegralType::cell); + kernels != _integrals.end()) + { + if (auto it = kernels->second.find(-1); it != kernels->second.end()) + { + std::vector& active_entities = it->second.second; + const int num_cells = topology.index_map(tdim)->size_local(); + active_entities.resize(num_cells); + std::iota(active_entities.begin(), active_entities.end(), 0); + } + } + + // Exterior facets. If there is a default integral, define it only + // on owned surface facets. + if (auto kernels = _integrals.find(IntegralType::exterior_facet); + kernels != _integrals.end()) + { + if (auto it = kernels->second.find(-1); it != kernels->second.end()) + { + std::vector& active_entities = it->second.second; + active_entities.clear(); + + // Get number of facets owned by this process + mesh.topology_mutable().create_connectivity(tdim - 1, tdim); + auto f_to_c = topology.connectivity(tdim - 1, tdim); + assert(topology.index_map(tdim - 1)); + std::set fwd_shared_facets; + + // Only need to consider shared facets when there are no ghost cells + if (topology.index_map(tdim)->num_ghosts() == 0) + { + fwd_shared_facets.insert( + topology.index_map(tdim - 1)->shared_indices().begin(), + topology.index_map(tdim - 1)->shared_indices().end()); + } + + const int num_facets = topology.index_map(tdim - 1)->size_local(); + for (int f = 0; f < num_facets; ++f) + { + if (f_to_c->num_links(f) == 1 + and fwd_shared_facets.find(f) == fwd_shared_facets.end()) + { + active_entities.push_back(f); + } + } + } + } + + // Interior facets. If there is a default integral, define it only on + // owned interior facets. + if (auto kernels = _integrals.find(IntegralType::interior_facet); + kernels != _integrals.end()) + { + if (auto it = kernels->second.find(-1); it != kernels->second.end()) + { + std::vector& active_entities = it->second.second; + + // Get number of facets owned by this process + mesh.topology_mutable().create_connectivity(tdim - 1, tdim); + assert(topology.index_map(tdim - 1)); + const int num_facets = topology.index_map(tdim - 1)->size_local(); + auto f_to_c = topology.connectivity(tdim - 1, tdim); + active_entities.clear(); + active_entities.reserve(num_facets); + for (int f = 0; f < num_facets; ++f) + { + if (f_to_c->num_links(f) == 2) + active_entities.push_back(f); + } + } + } + } // Function spaces (one for each argument) std::vector> _function_spaces; - // The mesh (needed for functionals when we don't have any spaces) + // Form coefficients + std::vector>> _coefficients; + + // Constants associated with the Form + std::vector>> _constants; + + // The mesh std::shared_ptr _mesh; -}; + + using kern + = std::function; + std::map>>> + _integrals; + + // True if permutation data needs to be passed into these integrals + bool _needs_permutation_data; + +}; // namespace fem } // namespace fem } // namespace dolfinx diff --git a/cpp/dolfinx/fem/FormCoefficients.h b/cpp/dolfinx/fem/FormCoefficients.h deleted file mode 100644 index 127e1d086d3..00000000000 --- a/cpp/dolfinx/fem/FormCoefficients.h +++ /dev/null @@ -1,132 +0,0 @@ -// Copyright (C) 2018 Chris Richardson -// -// This file is part of DOLFINX (https://www.fenicsproject.org) -// -// SPDX-License-Identifier: LGPL-3.0-or-later - -#pragma once - -#include -#include -#include -#include -#include -#include - -namespace dolfinx -{ - -namespace function -{ -template -class Function; -} - -namespace fem -{ -class FiniteElement; - -/// Storage for the coefficients of a Form consisting of Function and -/// the Element objects they are defined on. - -template -class FormCoefficients -{ -public: - /// Initialise the FormCoefficients, using tuples of - /// (original_coeff_position, name, Function). The Function pointer - /// may be a nullptr and assigned later. - FormCoefficients( - const std::vector< - std::tuple>>>& - coefficients) - { - for (const auto& c : coefficients) - { - _original_pos.push_back(std::get<0>(c)); - _names.push_back(std::get<1>(c)); - _coefficients.push_back(std::get<2>(c)); - } - } - - /// Get number of coefficients - int size() const { return _coefficients.size(); } - - /// Offset for each coefficient expansion array on a cell. Used to - /// pack data for multiple coefficients in a flat array. The last - /// entry is the size required to store all coefficients. - std::vector offsets() const - { - std::vector n{0}; - for (const auto& c : _coefficients) - { - if (!c) - throw std::runtime_error("Not all form coefficients have been set."); - n.push_back(n.back() + c->function_space()->element()->space_dimension()); - } - return n; - } - - /// Set coefficient with index i to be a Function - void set(int i, - const std::shared_ptr>& coefficient) - { - if (i >= (int)_coefficients.size()) - _coefficients.resize(i + 1); - _coefficients[i] = coefficient; - } - - /// Set coefficient with name to be a Function - void set(const std::string& name, - const std::shared_ptr>& coefficient) - { - const int i = get_index(name); - if (i >= (int)_coefficients.size()) - _coefficients.resize(i + 1); - _coefficients[i] = coefficient; - } - - /// Get the Function coefficient i - std::shared_ptr> get(int i) const - { - return _coefficients.at(i); - } - - /// Original position of coefficient in UFL form - /// @return The position of coefficient i in original ufl form - /// coefficients. - int original_position(int i) const { return _original_pos.at(i); } - - /// Get index from name of coefficient - /// @param[in] name Name of coefficient - /// @return Index of the coefficient - int get_index(const std::string& name) const - { - auto it = std::find(_names.begin(), _names.end(), name); - if (it == _names.end()) - throw std::runtime_error("Cannot find coefficient name:" + name); - return std::distance(_names.begin(), it); - } - - /// Get name from index of coefficient - /// @param[in] index Index of the coefficient - /// @return Name of the coefficient - std::string get_name(int index) const - { - if (index >= (int)_names.size()) - throw std::runtime_error("Invalid coefficient index"); - return _names[index]; - } - -private: - // Functions for the coefficients - std::vector>> _coefficients; - - // Copy of 'original positions' in UFL form - std::vector _original_pos; - - // Names of coefficients - std::vector _names; -}; -} // namespace fem -} // namespace dolfinx diff --git a/cpp/dolfinx/fem/FormIntegrals.h b/cpp/dolfinx/fem/FormIntegrals.h deleted file mode 100644 index 0b369e4ac13..00000000000 --- a/cpp/dolfinx/fem/FormIntegrals.h +++ /dev/null @@ -1,380 +0,0 @@ -// Copyright (C) 2018 Chris Richardson -// -// This file is part of DOLFINX (https://www.fenicsproject.org) -// -// SPDX-License-Identifier: LGPL-3.0-or-later - -#pragma once - -#include -#include -#include -#include -#include - -namespace dolfinx -{ -namespace mesh -{ -class Mesh; -} // namespace mesh - -namespace fem -{ - -/// Type of integral -enum class IntegralType : std::int8_t -{ - cell = 0, - exterior_facet = 1, - interior_facet = 2, - vertex = 3 -}; - -/// Integrals of a Form, including those defined over cells, interior -/// and exterior facets, and vertices. - -template -class FormIntegrals -{ -public: - /// Construct object from (index, tabulate function) pairs for - /// different integral types - /// @param[in] integrals For each integral type (domain index, - /// tabulate function) pairs for each integral. Domain index -1 means - /// for all entities. - /// @param[in] needs_permutation_data Pass true if an integral - /// requires mesh entity permutation data - FormIntegrals( - const std::map< - IntegralType, - std::vector>>>& integrals, - bool needs_permutation_data) - : _needs_permutation_data(needs_permutation_data) - { - for (auto& integral_type : integrals) - { - for (auto& integral : integral_type.second) - { - set_tabulate_tensor(integral_type.first, integral.first, - integral.second); - } - } - }; - - /// Get the function for 'tabulate_tensor' for integral i of given - /// type - /// @param[in] type Integral type - /// @param[in] i Integral number - /// @return Function to call for tabulate_tensor - const std::function& - get_tabulate_tensor(IntegralType type, int i) const - { - return _integrals.at(static_cast(type)).at(i).tabulate; - } - - /// @todo Should this be removed - /// - /// Set the function for 'tabulate_tensor' for integral i of - /// given type - /// @param[in] type Integral type - /// @param[in] i Integral number - /// @param[in] fn tabulate function - void set_tabulate_tensor( - IntegralType type, int i, - const std::function& fn) - { - std::vector& integrals - = _integrals.at(static_cast(type)); - - // Find insertion point - int pos = 0; - for (const auto& q : integrals) - { - if (q.id == i) - { - throw std::runtime_error("Integral with ID " + std::to_string(i) - + " already exists"); - } - else if (q.id > i) - break; - ++pos; - } - - // Insert new Integral - integrals.insert(integrals.begin() + pos, {fn, i, {}}); - } - - /// Get types of integrals in the form - /// @return Integrals types - std::set types() const - { - static const std::array types{IntegralType::cell, - IntegralType::exterior_facet, - IntegralType::interior_facet}; - std::set set; - for (auto type : types) - if (!_integrals.at(static_cast(type)).empty()) - set.insert(type); - return set; - } - - /// Number of integrals of given type - /// @param[in] type Integral type - /// @return Number of integrals - int num_integrals(IntegralType type) const - { - return _integrals.at(static_cast(type)).size(); - } - - /// Get the integer IDs of integrals of type t. The IDs correspond to - /// the domains which the integrals are defined for in the form, - /// except ID -1, which denotes the default integral. - /// @param[in] type Integral type - /// @return List of IDs for this integral - std::vector integral_ids(IntegralType type) const - { - std::vector ids; - for (const auto& integral : _integrals[static_cast(type)]) - ids.push_back(integral.id); - return ids; - } - - /// Get the list of active entities for the ith integral of type t. - /// Note, these are not retrieved by ID, but stored in order. The IDs - /// can be obtained with "FormIntegrals::integral_ids()". For cell - /// integrals, a list of cells. For facet integrals, a list of facets - /// etc. - /// @param[in] type Integral type - /// @param[in] i Integral number - /// @return List of active entities for this integral - const std::vector& integral_domains(IntegralType type, - int i) const - { - return _integrals.at(static_cast(type)).at(i).active_entities; - } - - /// Set the valid domains for the integrals of a given type from a - /// MeshTags "marker". Note the MeshTags is not stored, so if there - /// any changes to the integration domain this must be called again. - /// @param[in] type Integral type - /// @param[in] marker MeshTags mapping entities to integrals - void set_domains(IntegralType type, const mesh::MeshTags& marker) - { - std::vector& integrals - = _integrals.at(static_cast(type)); - if (integrals.size() == 0) - return; - - std::shared_ptr mesh = marker.mesh(); - const mesh::Topology& topology = mesh->topology(); - const int tdim = topology.dim(); - int dim = tdim; - if (type == IntegralType::exterior_facet - or type == IntegralType::interior_facet) - { - mesh->topology_mutable().create_connectivity(tdim - 1, tdim); - dim = tdim - 1; - } - else if (type == IntegralType::vertex) - dim = 0; - - if (dim != marker.dim()) - { - throw std::runtime_error("Invalid MeshTags dimension:" - + std::to_string(marker.dim())); - } - - // Create a reverse map - std::map id_to_integral; - for (std::size_t i = 0; i < integrals.size(); ++i) - { - if (integrals[i].id != -1) - { - integrals[i].active_entities.clear(); - id_to_integral.insert({integrals[i].id, i}); - } - } - - // Get mesh tag data - const std::vector& values = marker.values(); - const std::vector& tagged_entities = marker.indices(); - assert(topology.index_map(dim)); - const auto entity_end - = std::lower_bound(tagged_entities.begin(), tagged_entities.end(), - topology.index_map(dim)->size_local()); - - if (type == IntegralType::exterior_facet) - { - auto f_to_c = topology.connectivity(tdim - 1, tdim); - assert(f_to_c); - if (type == IntegralType::exterior_facet) - { - // Only need to consider shared facets when there are no ghost cells - assert(topology.index_map(tdim)); - std::set fwd_shared; - if (topology.index_map(tdim)->num_ghosts() == 0) - { - fwd_shared.insert( - topology.index_map(tdim - 1)->shared_indices().begin(), - topology.index_map(tdim - 1)->shared_indices().end()); - } - - for (auto f = tagged_entities.begin(); f != entity_end; ++f) - { - const std::size_t i = std::distance(tagged_entities.cbegin(), f); - - // All "owned" facets connected to one cell, that are not shared, - // should be external. - if (f_to_c->num_links(*f) == 1 - and fwd_shared.find(*f) == fwd_shared.end()) - { - if (auto it = id_to_integral.find(values[i]); - it != id_to_integral.end()) - { - integrals[it->second].active_entities.push_back(*f); - } - } - } - } - else if (type == IntegralType::interior_facet) - { - for (auto f = tagged_entities.begin(); f != entity_end; ++f) - { - if (f_to_c->num_links(*f) == 2) - { - const std::size_t i = std::distance(tagged_entities.cbegin(), f); - if (const auto it = id_to_integral.find(values[i]); - it != id_to_integral.end()) - { - integrals[it->second].active_entities.push_back(*f); - } - } - } - } - } - else - { - // For cell and vertex integrals use all markers (but not on ghost - // entities) - for (auto e = tagged_entities.begin(); e != entity_end; ++e) - { - const std::size_t i = std::distance(tagged_entities.cbegin(), e); - if (const auto it = id_to_integral.find(values[i]); - it != id_to_integral.end()) - { - integrals[it->second].active_entities.push_back(*e); - } - } - } - } - - /// If there exists a default integral of any type, set the list of - /// entities for those integrals from the mesh topology. For cell - /// integrals, this is all cells. For facet integrals, it is either - /// all interior or all exterior facets. - /// @param[in] mesh Mesh - void set_default_domains(const mesh::Mesh& mesh) - { - const mesh::Topology& topology = mesh.topology(); - const int tdim = topology.dim(); - std::vector& cell_integrals - = _integrals[static_cast(IntegralType::cell)]; - - // Cells. If there is a default integral, define it on all owned cells - if (cell_integrals.size() > 0 and cell_integrals[0].id == -1) - { - const int num_cells = topology.index_map(tdim)->size_local(); - cell_integrals[0].active_entities.resize(num_cells); - std::iota(cell_integrals[0].active_entities.begin(), - cell_integrals[0].active_entities.end(), 0); - } - - // Exterior facets. If there is a default integral, define it only on - // owned surface facets. - std::vector& exf_integrals - = _integrals[static_cast(IntegralType::exterior_facet)]; - if (exf_integrals.size() > 0 and exf_integrals[0].id == -1) - { - // If there is a default integral, define it only on surface facets - exf_integrals[0].active_entities.clear(); - - // Get number of facets owned by this process - mesh.topology_mutable().create_connectivity(tdim - 1, tdim); - auto f_to_c = topology.connectivity(tdim - 1, tdim); - assert(topology.index_map(tdim - 1)); - std::set fwd_shared_facets; - - // Only need to consider shared facets when there are no ghost cells - if (topology.index_map(tdim)->num_ghosts() == 0) - { - fwd_shared_facets.insert( - topology.index_map(tdim - 1)->shared_indices().begin(), - topology.index_map(tdim - 1)->shared_indices().end()); - } - - const int num_facets = topology.index_map(tdim - 1)->size_local(); - for (int f = 0; f < num_facets; ++f) - { - if (f_to_c->num_links(f) == 1 - and fwd_shared_facets.find(f) == fwd_shared_facets.end()) - exf_integrals[0].active_entities.push_back(f); - } - } - - // Interior facets. If there is a default integral, define it only on - // owned interior facets. - std::vector& inf_integrals - = _integrals[static_cast(IntegralType::interior_facet)]; - if (!inf_integrals.empty() and inf_integrals[0].id == -1) - { - // If there is a default integral, define it only on interior facets - inf_integrals[0].active_entities.clear(); - - // Get number of facets owned by this process - mesh.topology_mutable().create_connectivity(tdim - 1, tdim); - assert(topology.index_map(tdim - 1)); - const int num_facets = topology.index_map(tdim - 1)->size_local(); - auto f_to_c = topology.connectivity(tdim - 1, tdim); - inf_integrals[0].active_entities.reserve(num_facets); - for (int f = 0; f < num_facets; ++f) - { - if (f_to_c->num_links(f) == 2) - inf_integrals[0].active_entities.push_back(f); - } - } - } - - /// Get bool indicating whether permutation data needs to be passed - /// into these integrals - /// @return True if cell permutation data is required - bool needs_permutation_data() const { return _needs_permutation_data; } - -private: - // Collect together the function, id, and indices of entities to - // integrate on - struct Integral - { - std::function - tabulate; - int id; - std::vector active_entities; - }; - - // Array of vectors of integrals, arranged by type (see Type enum, and - // struct Integral above) - std::array, 4> _integrals; - - // A bool indicating whether permutation data needs to be passed into - // these integrals - bool _needs_permutation_data; -}; -} // namespace fem -} // namespace dolfinx diff --git a/cpp/dolfinx/fem/assemble_matrix_impl.h b/cpp/dolfinx/fem/assemble_matrix_impl.h index c8ce8a6d825..b8b44dbeebf 100644 --- a/cpp/dolfinx/fem/assemble_matrix_impl.h +++ b/cpp/dolfinx/fem/assemble_matrix_impl.h @@ -101,24 +101,23 @@ void assemble_matrix( = mesh->topology().connectivity(tdim, 0)->num_nodes(); // Get dofmap data - std::shared_ptr dofmap0 = a.function_space(0)->dofmap(); - std::shared_ptr dofmap1 = a.function_space(1)->dofmap(); + std::shared_ptr dofmap0 + = a.function_spaces().at(0)->dofmap(); + std::shared_ptr dofmap1 + = a.function_spaces().at(1)->dofmap(); assert(dofmap0); assert(dofmap1); const graph::AdjacencyList& dofs0 = dofmap0->list(); const graph::AdjacencyList& dofs1 = dofmap1->list(); // Prepare constants - if (!a.all_constants_set()) - throw std::runtime_error("Unset constant in Form"); const Eigen::Array constants = pack_constants(a); // Prepare coefficients const Eigen::Array coeffs = pack_coefficients(a); - const FormIntegrals& integrals = a.integrals(); - const bool needs_permutation_data = integrals.needs_permutation_data(); + const bool needs_permutation_data = a.needs_permutation_data(); if (needs_permutation_data) mesh->topology_mutable().create_entity_permutations(); const Eigen::Array& cell_info @@ -126,18 +125,18 @@ void assemble_matrix( ? mesh->topology().get_cell_permutation_info() : Eigen::Array(num_cells); - for (int i = 0; i < integrals.num_integrals(IntegralType::cell); ++i) + for (int i : a.integral_ids(IntegralType::cell)) { - const auto& fn = integrals.get_tabulate_tensor(IntegralType::cell, i); + const auto& fn = a.kernel(IntegralType::cell, i); const std::vector& active_cells - = integrals.integral_domains(IntegralType::cell, i); + = a.domains(IntegralType::cell, i); impl::assemble_cells(mat_set_values, mesh->geometry(), active_cells, dofs0, dofs1, bc0, bc1, fn, coeffs, constants, cell_info); } - if (integrals.num_integrals(IntegralType::exterior_facet) > 0 - or integrals.num_integrals(IntegralType::interior_facet) > 0) + if (a.num_integrals(IntegralType::exterior_facet) > 0 + or a.num_integrals(IntegralType::interior_facet) > 0) { // FIXME: cleanup these calls? Some of the happen internally again. mesh->topology_mutable().create_entities(tdim - 1); @@ -151,26 +150,22 @@ void assemble_matrix( : Eigen::Array( facets_per_cell, num_cells); - for (int i = 0; i < integrals.num_integrals(IntegralType::exterior_facet); - ++i) + for (int i : a.integral_ids(IntegralType::exterior_facet)) { - const auto& fn - = integrals.get_tabulate_tensor(IntegralType::exterior_facet, i); + const auto& fn = a.kernel(IntegralType::exterior_facet, i); const std::vector& active_facets - = integrals.integral_domains(IntegralType::exterior_facet, i); + = a.domains(IntegralType::exterior_facet, i); impl::assemble_exterior_facets(mat_set_values, *mesh, active_facets, dofs0, dofs1, bc0, bc1, fn, coeffs, constants, cell_info, perms); } - const std::vector c_offsets = a.coefficients().offsets(); - for (int i = 0; i < integrals.num_integrals(IntegralType::interior_facet); - ++i) + const std::vector c_offsets = a.coefficient_offsets(); + for (int i : a.integral_ids(IntegralType::interior_facet)) { - const auto& fn - = integrals.get_tabulate_tensor(IntegralType::interior_facet, i); + const auto& fn = a.kernel(IntegralType::interior_facet, i); const std::vector& active_facets - = integrals.integral_domains(IntegralType::interior_facet, i); + = a.domains(IntegralType::interior_facet, i); impl::assemble_interior_facets( mat_set_values, *mesh, active_facets, *dofmap0, *dofmap1, bc0, bc1, fn, coeffs, c_offsets, constants, cell_info, perms); diff --git a/cpp/dolfinx/fem/assemble_scalar_impl.h b/cpp/dolfinx/fem/assemble_scalar_impl.h index f42bb375f60..82a13f19c65 100644 --- a/cpp/dolfinx/fem/assemble_scalar_impl.h +++ b/cpp/dolfinx/fem/assemble_scalar_impl.h @@ -73,15 +73,14 @@ T assemble_scalar(const fem::Form& M) = mesh->topology().connectivity(tdim, 0)->num_nodes(); // Prepare constants - if (!M.all_constants_set()) - throw std::runtime_error("Unset constant in Form"); - auto constants = M.constants(); + const std::vector>>& constants + = M.constants(); std::vector constant_values; for (auto const& constant : constants) { // Get underlying data array of this Constant - const std::vector& array = constant.second->value; + const std::vector& array = constant->value; constant_values.insert(constant_values.end(), array.data(), array.data() + array.size()); } @@ -90,8 +89,7 @@ T assemble_scalar(const fem::Form& M) const Eigen::Array coeffs = pack_coefficients(M); - const FormIntegrals& integrals = M.integrals(); - const bool needs_permutation_data = integrals.needs_permutation_data(); + const bool needs_permutation_data = M.needs_permutation_data(); if (needs_permutation_data) mesh->topology_mutable().create_entity_permutations(); const Eigen::Array& cell_info @@ -100,17 +98,17 @@ T assemble_scalar(const fem::Form& M) : Eigen::Array(num_cells); T value(0); - for (int i = 0; i < integrals.num_integrals(IntegralType::cell); ++i) + for (int i : M.integral_ids(IntegralType::cell)) { - const auto& fn = integrals.get_tabulate_tensor(IntegralType::cell, i); + const auto& fn = M.kernel(IntegralType::cell, i); const std::vector& active_cells - = integrals.integral_domains(IntegralType::cell, i); + = M.domains(IntegralType::cell, i); value += fem::impl::assemble_cells(mesh->geometry(), active_cells, fn, coeffs, constant_values, cell_info); } - if (integrals.num_integrals(IntegralType::exterior_facet) > 0 - or integrals.num_integrals(IntegralType::interior_facet) > 0) + if (M.num_integrals(IntegralType::exterior_facet) > 0 + or M.num_integrals(IntegralType::interior_facet) > 0) { // FIXME: cleanup these calls? Some of these happen internally again. mesh->topology_mutable().create_entities(tdim - 1); @@ -124,25 +122,21 @@ T assemble_scalar(const fem::Form& M) : Eigen::Array( facets_per_cell, num_cells); - for (int i = 0; i < integrals.num_integrals(IntegralType::exterior_facet); - ++i) + for (int i : M.integral_ids(IntegralType::exterior_facet)) { - const auto& fn - = integrals.get_tabulate_tensor(IntegralType::exterior_facet, i); + const auto& fn = M.kernel(IntegralType::exterior_facet, i); const std::vector& active_facets - = integrals.integral_domains(IntegralType::exterior_facet, i); + = M.domains(IntegralType::exterior_facet, i); value += fem::impl::assemble_exterior_facets( *mesh, active_facets, fn, coeffs, constant_values, cell_info, perms); } - const std::vector c_offsets = M.coefficients().offsets(); - for (int i = 0; i < integrals.num_integrals(IntegralType::interior_facet); - ++i) + const std::vector c_offsets = M.coefficient_offsets(); + for (int i : M.integral_ids(IntegralType::interior_facet)) { - const auto& fn - = integrals.get_tabulate_tensor(IntegralType::interior_facet, i); + const auto& fn = M.kernel(IntegralType::interior_facet, i); const std::vector& active_facets - = integrals.integral_domains(IntegralType::interior_facet, i); + = M.domains(IntegralType::interior_facet, i); value += fem::impl::assemble_interior_facets( *mesh, active_facets, fn, coeffs, c_offsets, constant_values, cell_info, perms); diff --git a/cpp/dolfinx/fem/assemble_vector_impl.h b/cpp/dolfinx/fem/assemble_vector_impl.h index 0c1e1f27a4b..2b75a9f5114 100644 --- a/cpp/dolfinx/fem/assemble_vector_impl.h +++ b/cpp/dolfinx/fem/assemble_vector_impl.h @@ -158,10 +158,10 @@ void _lift_bc_cells( mesh->topology_mutable().create_entity_permutations(); // Get dofmap for columns and rows of a - assert(a.function_space(0)); - assert(a.function_space(1)); - std::shared_ptr dofmap0 = a.function_space(0)->dofmap(); - std::shared_ptr dofmap1 = a.function_space(1)->dofmap(); + assert(a.function_spaces().at(0)); + assert(a.function_spaces().at(1)); + std::shared_ptr dofmap0 = a.function_spaces()[0]->dofmap(); + std::shared_ptr dofmap1 = a.function_spaces()[1]->dofmap(); assert(dofmap0); assert(dofmap1); @@ -171,7 +171,7 @@ void _lift_bc_cells( const std::function& fn - = a.integrals().get_tabulate_tensor(IntegralType::cell, 0); + = a.kernel(IntegralType::cell, -1); // Prepare cell geometry const int gdim = mesh->geometry().dim(); @@ -190,8 +190,6 @@ void _lift_bc_cells( Eigen::Matrix be; // Prepare constants - if (!a.all_constants_set()) - throw std::runtime_error("Unset constant in Form"); const Eigen::Array constant_values = pack_constants(a); const Eigen::Array& cell_info @@ -281,10 +279,12 @@ void _lift_bc_exterior_facets( mesh->topology_mutable().create_entity_permutations(); // Get dofmap for columns and rows of a - assert(a.function_space(0)); - assert(a.function_space(1)); - std::shared_ptr dofmap0 = a.function_space(0)->dofmap(); - std::shared_ptr dofmap1 = a.function_space(1)->dofmap(); + assert(a.function_spaces().at(0)); + assert(a.function_spaces().at(1)); + std::shared_ptr dofmap0 + = a.function_spaces().at(0)->dofmap(); + std::shared_ptr dofmap1 + = a.function_spaces().at(1)->dofmap(); assert(dofmap0); assert(dofmap1); @@ -294,7 +294,7 @@ void _lift_bc_exterior_facets( const std::function& fn - = a.integrals().get_tabulate_tensor(IntegralType::exterior_facet, 0); + = a.kernel(IntegralType::exterior_facet, -1); // Prepare cell geometry const graph::AdjacencyList& x_dofmap @@ -311,8 +311,6 @@ void _lift_bc_exterior_facets( Eigen::Matrix be; // Prepare constants - if (!a.all_constants_set()) - throw std::runtime_error("Unset constant in Form"); const Eigen::Array constant_values = pack_constants(a); // Iterate over owned facets @@ -425,22 +423,20 @@ void assemble_vector(Eigen::Ref> b, = mesh->topology().connectivity(tdim, 0)->num_nodes(); // Get dofmap data - assert(L.function_space(0)); - std::shared_ptr dofmap = L.function_space(0)->dofmap(); + assert(L.function_spaces().at(0)); + std::shared_ptr dofmap + = L.function_spaces().at(0)->dofmap(); assert(dofmap); const graph::AdjacencyList& dofs = dofmap->list(); // Prepare constants - if (!L.all_constants_set()) - throw std::runtime_error("Unset constant in Form"); const Eigen::Array constant_values = pack_constants(L); // Prepare coefficients const Eigen::Array coeffs = pack_coefficients(L); - const FormIntegrals& integrals = L.integrals(); - const bool needs_permutation_data = integrals.needs_permutation_data(); + const bool needs_permutation_data = L.needs_permutation_data(); if (needs_permutation_data) mesh->topology_mutable().create_entity_permutations(); const Eigen::Array& cell_info @@ -448,17 +444,17 @@ void assemble_vector(Eigen::Ref> b, ? mesh->topology().get_cell_permutation_info() : Eigen::Array(num_cells); - for (int i = 0; i < integrals.num_integrals(IntegralType::cell); ++i) + for (int i : L.integral_ids(IntegralType::cell)) { - const auto& fn = integrals.get_tabulate_tensor(IntegralType::cell, i); + const auto& fn = L.kernel(IntegralType::cell, i); const std::vector& active_cells - = integrals.integral_domains(IntegralType::cell, i); + = L.domains(IntegralType::cell, i); fem::impl::assemble_cells(b, mesh->geometry(), active_cells, dofs, fn, coeffs, constant_values, cell_info); } - if (integrals.num_integrals(IntegralType::exterior_facet) > 0 - or integrals.num_integrals(IntegralType::interior_facet) > 0) + if (L.num_integrals(IntegralType::exterior_facet) > 0 + or L.num_integrals(IntegralType::interior_facet) > 0) { // FIXME: cleanup these calls? Some of the happen internally again. mesh->topology_mutable().create_entities(tdim - 1); @@ -471,26 +467,22 @@ void assemble_vector(Eigen::Ref> b, ? mesh->topology().get_facet_permutations() : Eigen::Array( facets_per_cell, num_cells); - for (int i = 0; i < integrals.num_integrals(IntegralType::exterior_facet); - ++i) + for (int i : L.integral_ids(IntegralType::exterior_facet)) { - const auto& fn - = integrals.get_tabulate_tensor(IntegralType::exterior_facet, i); + const auto& fn = L.kernel(IntegralType::exterior_facet, i); const std::vector& active_facets - = integrals.integral_domains(IntegralType::exterior_facet, i); + = L.domains(IntegralType::exterior_facet, i); fem::impl::assemble_exterior_facets(b, *mesh, active_facets, dofs, fn, coeffs, constant_values, cell_info, perms); } - const std::vector c_offsets = L.coefficients().offsets(); - for (int i = 0; i < integrals.num_integrals(IntegralType::interior_facet); - ++i) + const std::vector c_offsets = L.coefficient_offsets(); + for (int i : L.integral_ids(IntegralType::interior_facet)) { - const auto& fn - = integrals.get_tabulate_tensor(IntegralType::interior_facet, i); + const auto& fn = L.kernel(IntegralType::interior_facet, i); const std::vector& active_facets - = integrals.integral_domains(IntegralType::interior_facet, i); + = L.domains(IntegralType::interior_facet, i); fem::impl::assemble_interior_facets(b, *mesh, active_facets, *dofmap, fn, coeffs, c_offsets, constant_values, cell_info, perms); @@ -744,7 +736,7 @@ void apply_lifting( Eigen::Matrix bc_values1; if (a[j] and !bcs1[j].empty()) { - auto V1 = a[j]->function_space(1); + auto V1 = a[j]->function_spaces()[1]; assert(V1); auto map1 = V1->dofmap()->index_map; assert(map1); @@ -774,9 +766,9 @@ void lift_bc( const std::vector& bc_markers1, double scale) { const Eigen::Matrix x0(0); - if (a.integrals().num_integrals(fem::IntegralType::cell) > 0) + if (a.num_integrals(fem::IntegralType::cell) > 0) _lift_bc_cells(b, a, bc_values1, bc_markers1, x0, scale); - if (a.integrals().num_integrals(fem::IntegralType::exterior_facet) > 0) + if (a.num_integrals(fem::IntegralType::exterior_facet) > 0) _lift_bc_exterior_facets(b, a, bc_values1, bc_markers1, x0, scale); } //----------------------------------------------------------------------------- @@ -788,9 +780,9 @@ void lift_bc( const Eigen::Ref>& x0, double scale) { - if (a.integrals().num_integrals(fem::IntegralType::cell) > 0) + if (a.num_integrals(fem::IntegralType::cell) > 0) _lift_bc_cells(b, a, bc_values1, bc_markers1, x0, scale); - if (a.integrals().num_integrals(fem::IntegralType::exterior_facet) > 0) + if (a.num_integrals(fem::IntegralType::exterior_facet) > 0) _lift_bc_exterior_facets(b, a, bc_values1, bc_markers1, x0, scale); } //----------------------------------------------------------------------------- diff --git a/cpp/dolfinx/fem/assembler.h b/cpp/dolfinx/fem/assembler.h index 3d0d24ea5a1..bd029c45428 100644 --- a/cpp/dolfinx/fem/assembler.h +++ b/cpp/dolfinx/fem/assembler.h @@ -100,25 +100,27 @@ void assemble_matrix( { // Index maps for dof ranges - auto map0 = a.function_space(0)->dofmap()->index_map; - auto map1 = a.function_space(1)->dofmap()->index_map; + auto map0 = a.function_spaces().at(0)->dofmap()->index_map; + auto map1 = a.function_spaces().at(1)->dofmap()->index_map; // Build dof markers std::vector dof_marker0, dof_marker1; + assert(map0); std::int32_t dim0 = map0->block_size() * (map0->size_local() + map0->num_ghosts()); + assert(map1); std::int32_t dim1 = map1->block_size() * (map1->size_local() + map1->num_ghosts()); for (std::size_t k = 0; k < bcs.size(); ++k) { assert(bcs[k]); assert(bcs[k]->function_space()); - if (a.function_space(0)->contains(*bcs[k]->function_space())) + if (a.function_spaces().at(0)->contains(*bcs[k]->function_space())) { dof_marker0.resize(dim0, false); bcs[k]->mark_dofs(dof_marker0); } - if (a.function_space(1)->contains(*bcs[k]->function_space())) + if (a.function_spaces().at(1)->contains(*bcs[k]->function_space())) { dof_marker1.resize(dim1, false); bcs[k]->mark_dofs(dof_marker1); @@ -265,7 +267,7 @@ bcs_rows(const std::vector*>& L, L.size()); for (std::size_t i = 0; i < L.size(); ++i) for (const std::shared_ptr>& bc : bcs) - if (L[i]->function_space(0)->contains(*bc->function_space())) + if (L[i]->function_spaces()[0]->contains(*bc->function_space())) bcs0[i].push_back(bc); return bcs0; } @@ -300,7 +302,7 @@ bcs_cols(const std::vector>>>& a, // FIXME: handle case where a[i][j] is null if (a[i][j]) { - if (a[i][j]->function_space(1)->contains(*bc->function_space())) + if (a[i][j]->function_spaces()[1]->contains(*bc->function_space())) bcs1[i][j].push_back(bc); } } diff --git a/cpp/dolfinx/fem/petsc.cpp b/cpp/dolfinx/fem/petsc.cpp index 772d3041f94..f6621099f48 100644 --- a/cpp/dolfinx/fem/petsc.cpp +++ b/cpp/dolfinx/fem/petsc.cpp @@ -62,16 +62,16 @@ la::PETScMatrix fem::create_matrix_block( assert(patterns[row].back()); auto& sp = patterns[row].back(); assert(sp); - const FormIntegrals& integrals = a(row, col)->integrals(); - if (integrals.num_integrals(IntegralType::cell) > 0) + const fem::Form& a_ = *a(row, col); + if (a_.num_integrals(IntegralType::cell) > 0) SparsityPatternBuilder::cells(*sp, mesh->topology(), dofmaps); - if (integrals.num_integrals(IntegralType::interior_facet) > 0) + if (a_.num_integrals(IntegralType::interior_facet) > 0) { mesh->topology_mutable().create_entities(tdim - 1); SparsityPatternBuilder::interior_facets(*sp, mesh->topology(), dofmaps); } - if (integrals.num_integrals(IntegralType::exterior_facet) > 0) + if (a_.num_integrals(IntegralType::exterior_facet) > 0) { mesh->topology_mutable().create_entities(tdim - 1); SparsityPatternBuilder::exterior_facets(*sp, mesh->topology(), diff --git a/cpp/dolfinx/fem/utils.cpp b/cpp/dolfinx/fem/utils.cpp index 2a433d5221f..67596ab563e 100644 --- a/cpp/dolfinx/fem/utils.cpp +++ b/cpp/dolfinx/fem/utils.cpp @@ -206,6 +206,24 @@ fem::DofMap fem::create_dofmap(MPI_Comm comm, const ufc_dofmap& ufc_dofmap, return DofMap(element_dof_layout, index_map, std::move(dofmap)); } //----------------------------------------------------------------------------- +std::vector fem::get_coefficient_names(const ufc_form& ufc_form) +{ + std::vector coefficients; + const char** names = ufc_form.coefficient_name_map(); + for (int i = 0; i < ufc_form.num_coefficients; ++i) + coefficients.push_back(names[i]); + return coefficients; +} +//----------------------------------------------------------------------------- +std::vector fem::get_constant_names(const ufc_form& ufc_form) +{ + std::vector constants; + const char** names = ufc_form.constant_name_map(); + for (int i = 0; i < ufc_form.num_constants; ++i) + constants.push_back(names[i]); + return constants; +} +//----------------------------------------------------------------------------- fem::CoordinateElement fem::create_coordinate_map(const ufc_coordinate_mapping& ufc_cmap) { diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 022c8664549..2251d14885d 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -67,9 +67,16 @@ extract_function_spaces( std::array, 2>>( a.cols())); for (int i = 0; i < a.rows(); ++i) + { for (int j = 0; j < a.cols(); ++j) + { if (a(i, j)) - spaces[i][j] = {a(i, j)->function_space(0), a(i, j)->function_space(1)}; + { + spaces[i][j] + = {a(i, j)->function_spaces()[0], a(i, j)->function_spaces()[1]}; + } + } + } return spaces; } @@ -88,12 +95,12 @@ la::SparsityPattern create_sparsity_pattern(const Form& a) } // Get dof maps and mesh - std::array dofmaps{a.function_space(0)->dofmap().get(), - a.function_space(1)->dofmap().get()}; + std::array dofmaps{a.function_spaces().at(0)->dofmap().get(), + a.function_spaces().at(1)->dofmap().get()}; std::shared_ptr mesh = a.mesh(); assert(mesh); - const std::set types = a.integrals().types(); + const std::set types = a.integral_types(); if (types.find(IntegralType::interior_facet) != types.end() or types.find(IntegralType::exterior_facet) != types.end()) { @@ -127,50 +134,48 @@ ElementDofLayout create_element_dof_layout(const ufc_dofmap& dofmap, DofMap create_dofmap(MPI_Comm comm, const ufc_dofmap& dofmap, mesh::Topology& topology); -/// Extract coefficients from a UFC form -template -std::vector< - std::tuple>>> -get_coeffs_from_ufc_form(const ufc_form& ufc_form) -{ - std::vector< - std::tuple>>> - coeffs; - const char** names = ufc_form.coefficient_name_map(); - for (int i = 0; i < ufc_form.num_coefficients; ++i) - { - coeffs.emplace_back(ufc_form.original_coefficient_position(i), names[i], - nullptr); - } - return coeffs; -} +/// Get the name of each coefficient in a UFC form +/// @param[in] ufc_form The UFC form +/// return The name of each coefficient +std::vector get_coefficient_names(const ufc_form& ufc_form); -/// Extract coefficients from a UFC form -template -std::vector< - std::pair>>> -get_constants_from_ufc_form(const ufc_form& ufc_form) -{ - std::vector< - std::pair>>> - constants; - const char** names = ufc_form.constant_name_map(); - for (int i = 0; i < ufc_form.num_constants; ++i) - constants.emplace_back(names[i], nullptr); - return constants; -} +/// Get the name of each constant in a UFC form +/// @param[in] ufc_form The UFC form +/// return The name of each constant +std::vector get_constant_names(const ufc_form& ufc_form); /// Create a Form from UFC input /// @param[in] ufc_form The UFC form /// @param[in] spaces Vector of function spaces +/// @param[in] coefficients Coefficient fields in the form +/// @param[in] constants Spatial constants in the form +/// @param[in] subdomains Subdomain markers +/// @param[in] mesh The mesh of the domain template Form create_form( const ufc_form& ufc_form, - const std::vector>& spaces) + const std::vector>& spaces, + const std::vector>>& + coefficients, + const std::vector>>& constants, + const std::map*>& subdomains, + const std::shared_ptr& mesh = nullptr) { - assert(ufc_form.rank == (int)spaces.size()); + if (ufc_form.rank != (int)spaces.size()) + throw std::runtime_error("Wrong number of argument spaces for Form."); + if (ufc_form.num_coefficients != (int)coefficients.size()) + { + throw std::runtime_error( + "Mismatch between number of expected and provided Form coefficients."); + } + if (ufc_form.num_constants != (int)constants.size()) + { + throw std::runtime_error( + "Mismatch between number of expected and provided Form constants."); + } // Check argument function spaces +#ifdef DEBUG for (std::size_t i = 0; i < spaces.size(); ++i) { assert(spaces[i]->element()); @@ -184,16 +189,19 @@ Form create_form( "Cannot create form. Wrong type of function space for argument."); } } +#endif // Get list of integral IDs, and load tabulate tensor into memory for each - bool needs_permutation_data = false; - std::map>>> + using kern = std::function; + std::map>, + const mesh::MeshTags*>> integral_data; + bool needs_permutation_data = false; + + // Attach cell kernels std::vector cell_integral_ids(ufc_form.num_cell_integrals); ufc_form.get_cell_integral_ids(cell_integral_ids.data()); for (int id : cell_integral_ids) @@ -202,13 +210,21 @@ Form create_form( assert(integral); if (integral->needs_permutation_data) needs_permutation_data = true; - integral_data[IntegralType::cell].emplace_back(id, - integral->tabulate_tensor); + integral_data[IntegralType::cell].first.emplace_back( + id, integral->tabulate_tensor); std::free(integral); } - // FIXME: Can this be handled better? - // FIXME: Handle forms with no space + // Attach cell subdomain data + if (auto it = subdomains.find(IntegralType::cell); + it != subdomains.end() and !cell_integral_ids.empty()) + { + integral_data[IntegralType::cell].second = it->second; + } + + // FIXME: Can facets be handled better? + + // Create facets, if required if (ufc_form.num_exterior_facet_integrals > 0 or ufc_form.num_interior_facet_integrals > 0) { @@ -220,6 +236,7 @@ Form create_form( } } + // Attach exterior facet kernels std::vector exterior_facet_integral_ids( ufc_form.num_exterior_facet_integrals); ufc_form.get_exterior_facet_integral_ids(exterior_facet_integral_ids.data()); @@ -229,11 +246,19 @@ Form create_form( assert(integral); if (integral->needs_permutation_data) needs_permutation_data = true; - integral_data[IntegralType::exterior_facet].emplace_back( + integral_data[IntegralType::exterior_facet].first.emplace_back( id, integral->tabulate_tensor); std::free(integral); } + // Attach exterior facet subdomain data + if (auto it = subdomains.find(IntegralType::exterior_facet); + it != subdomains.end() and !exterior_facet_integral_ids.empty()) + { + integral_data[IntegralType::exterior_facet].second = it->second; + } + + // Attach interior facet kernels std::vector interior_facet_integral_ids( ufc_form.num_interior_facet_integrals); ufc_form.get_interior_facet_integral_ids(interior_facet_integral_ids.data()); @@ -243,40 +268,104 @@ Form create_form( assert(integral); if (integral->needs_permutation_data) needs_permutation_data = true; - integral_data[IntegralType::interior_facet].emplace_back( + integral_data[IntegralType::interior_facet].first.emplace_back( id, integral->tabulate_tensor); std::free(integral); } - // Not currently working + // Attach interior facet subdomain data + if (auto it = subdomains.find(IntegralType::interior_facet); + it != subdomains.end() and !interior_facet_integral_ids.empty()) + { + integral_data[IntegralType::interior_facet].second = it->second; + } + + // Vertex integrals: not currently working std::vector vertex_integral_ids(ufc_form.num_vertex_integrals); ufc_form.get_vertex_integral_ids(vertex_integral_ids.data()); - if (vertex_integral_ids.size() > 0) + if (!vertex_integral_ids.empty()) { throw std::runtime_error( "Vertex integrals not supported. Under development."); } - return fem::Form( - spaces, FormIntegrals(integral_data, needs_permutation_data), - FormCoefficients(fem::get_coeffs_from_ufc_form(ufc_form)), - fem::get_constants_from_ufc_form(ufc_form)); + return fem::Form(spaces, integral_data, coefficients, constants, + needs_permutation_data, mesh); +} + +/// Create a Form from UFC input +/// @param[in] ufc_form The UFC form +/// @param[in] spaces The function spaces for the Form arguments +/// @param[in] coefficients Coefficient fields in the form (by name) +/// @param[in] constants Spatial constants in the form (by name) +/// @param[in] subdomains Subdomain makers +/// @param[in] mesh The mesh of the domain. This is required if the form +/// has no arguments, e.g. a functional. +/// @return A Form +template +Form create_form( + const ufc_form& ufc_form, + const std::vector>& spaces, + const std::map>>& + coefficients, + const std::map>>& + constants, + const std::map*>& subdomains, + const std::shared_ptr& mesh = nullptr) +{ + // Place coefficients in appropriate order + std::vector>> coeff_map; + for (const std::string& name : get_coefficient_names(ufc_form)) + { + if (auto it = coefficients.find(name); it != coefficients.end()) + coeff_map.push_back(it->second); + else + { + throw std::runtime_error("Form coefficient \"" + name + + "\" not provided."); + } + } + + // Place constants in appropriate order + std::vector>> const_map; + for (const std::string& name : get_constant_names(ufc_form)) + { + if (auto it = constants.find(name); it != constants.end()) + const_map.push_back(it->second); + else + { + throw std::runtime_error("Form constant \"" + name + "\" not provided."); + } + } + + return create_form(ufc_form, spaces, coeff_map, const_map, subdomains, mesh); } -/// Create a form from a form_create function returning a pointer to a -/// ufc_form, taking care of memory allocation +/// Create a Form using a factory function that returns a pointer to a +/// ufc_form. /// @param[in] fptr pointer to a function returning a pointer to -/// ufc_form -/// @param[in] spaces function spaces -/// @return Form +/// ufc_form +/// @param[in] spaces The function spaces for the Form arguments +/// @param[in] coefficients Coefficient fields in the form (by name) +/// @param[in] constants Spatial constants in the form (by name) +/// @param[in] subdomains Subdomain markers +/// @param[in] mesh The mesh of the domain. This is required if the form +/// has no arguments, e.g. a functional. +/// @return A Form template std::shared_ptr> create_form( ufc_form* (*fptr)(), - const std::vector>& spaces) + const std::vector>& spaces, + const std::map>>& + coefficients, + const std::map>>& + constants, + const std::map*>& subdomains, + const std::shared_ptr& mesh = nullptr) { ufc_form* form = fptr(); - auto L = std::make_shared>( - dolfinx::fem::create_form(*form, spaces)); + auto L = std::make_shared>(dolfinx::fem::create_form( + *form, spaces, coefficients, constants, subdomains, mesh)); std::free(form); return L; } @@ -314,14 +403,15 @@ Eigen::Array pack_coefficients(const fem::Form& form) { // Get form coefficient offsets amd dofmaps - const fem::FormCoefficients& coefficients = form.coefficients(); - const std::vector& offsets = coefficients.offsets(); + const std::vector>> coefficients + = form.coefficients(); + const std::vector offsets = form.coefficient_offsets(); std::vector dofmaps(coefficients.size()); std::vector>> v; - for (int i = 0; i < coefficients.size(); ++i) + for (std::size_t i = 0; i < coefficients.size(); ++i) { - dofmaps[i] = coefficients.get(i)->function_space()->dofmap().get(); - v.emplace_back(coefficients.get(i)->x()->array()); + dofmaps[i] = coefficients[i]->function_space()->dofmap().get(); + v.emplace_back(coefficients[i]->x()->array()); } // Get mesh @@ -361,7 +451,7 @@ Eigen::Array pack_constants(const fem::Form& form) std::vector constant_values; for (auto& constant : form.constants()) { - const std::vector& array = constant.second->value; + const std::vector& array = constant->value; constant_values.insert(constant_values.end(), array.begin(), array.end()); } diff --git a/python/demo/mixed-elasticity-sc/static-condensation-elasticity.py b/python/demo/mixed-elasticity-sc/static-condensation-elasticity.py index d71034f7618..78861a64361 100644 --- a/python/demo/mixed-elasticity-sc/static-condensation-elasticity.py +++ b/python/demo/mixed-elasticity-sc/static-condensation-elasticity.py @@ -140,9 +140,9 @@ def tabulate_condensed_tensor_A(A_, w_, c_, coords_, entity_local_index, permuta A[:, :] = - A10 @ numpy.linalg.solve(A00, A01) -# Prepare an empty Form and set the condensed tabulation kernel -a_cond = dolfinx.cpp.fem.Form([U._cpp_object, U._cpp_object], False) -a_cond.set_tabulate_tensor(dolfinx.fem.IntegralType.cell, -1, tabulate_condensed_tensor_A.address) +# Prepare a Form with a condensed tabulation kernel +integrals = {dolfinx.fem.IntegralType.cell: ([(-1, tabulate_condensed_tensor_A.address)], None)} +a_cond = dolfinx.cpp.fem.Form([U._cpp_object, U._cpp_object], integrals, [], [], False, None) A_cond = dolfinx.fem.assemble_matrix(a_cond, [bc]) A_cond.assemble() diff --git a/python/dolfinx/fem/assemble.py b/python/dolfinx/fem/assemble.py index acbf8a193a5..738e8b7ce29 100644 --- a/python/dolfinx/fem/assemble.py +++ b/python/dolfinx/fem/assemble.py @@ -288,7 +288,7 @@ def fn(form): return rows, cols -@ assemble_matrix_block.register(PETSc.Mat) +@assemble_matrix_block.register(PETSc.Mat) def _(A: PETSc.Mat, a: typing.List[typing.List[typing.Union[Form, cpp.fem.Form]]], bcs: typing.List[DirichletBC] = [], diff --git a/python/dolfinx/fem/form.py b/python/dolfinx/fem/form.py index 93fceed5c80..7db08783687 100644 --- a/python/dolfinx/fem/form.py +++ b/python/dolfinx/fem/form.py @@ -37,6 +37,8 @@ def __init__(self, form: ufl.Form, form_compiler_parameters: dict = {}, jit_para self._subdomains, = list(sd.values()) # Assuming single domain domain, = list(sd.keys()) # Assuming single domain mesh = domain.ufl_cargo() + if mesh is None: + raise RuntimeError("Expecting to find a Mesh in the form.") # Compile UFL form with JIT ufc_form = jit.ffcx_jit( @@ -50,43 +52,21 @@ def __init__(self, form: ufl.Form, form_compiler_parameters: dict = {}, jit_para func.ufl_function_space()._cpp_object for func in form.arguments() ] - # Prepare dolfinx.cpp.fem.Form and hold it as a member - ffi = cffi.FFI() - self._cpp_object = cpp.fem.create_form(ffi.cast("uintptr_t", ufc_form), function_spaces) - - # Need to fill the form with coefficients data - # For every coefficient in form take its C++ object + # Prepare coefficients data. For every coefficient in form take + # its C++ object. original_coefficients = form.coefficients() - for i in range(self._cpp_object.num_coefficients()): - j = self._cpp_object.original_coefficient_position(i) - self._cpp_object.set_coefficient(i, original_coefficients[j]._cpp_object) - - # Constants are set based on their position in original form - original_constants = [c._cpp_object for c in form.constants()] + coeffs = [original_coefficients[ufc_form.original_coefficient_position( + i)]._cpp_object for i in range(ufc_form.num_coefficients)] - self._cpp_object.set_constants(original_constants) + # Create dictionary of of subdomain markers (possible None for + # some dimensions + subdomains = {cpp.fem.IntegralType.cell: self._subdomains.get("cell"), + cpp.fem.IntegralType.exterior_facet: self._subdomains.get("exterior_facet"), + cpp.fem.IntegralType.interior_facet: self._subdomains.get("interior_facet"), + cpp.fem.IntegralType.vertex: self._subdomains.get("vertex")} - if mesh is None: - raise RuntimeError("Expecting to find a Mesh in the form.") - - # Attach mesh (because function spaces and coefficients may be - # empty lists) - if not function_spaces: - self._cpp_object.set_mesh(mesh) - - # Attach subdomains to C++ Form if we have them - subdomains = self._subdomains.get("cell") - if subdomains: - self._cpp_object.integrals.set_domains(cpp.fem.IntegralType.cell, subdomains) - - subdomains = self._subdomains.get("exterior_facet") - if subdomains: - self._cpp_object.integrals.set_domains(cpp.fem.IntegralType.exterior_facet, subdomains) - - subdomains = self._subdomains.get("interior_facet") - if subdomains: - self._cpp_object.integrals.set_domains(cpp.fem.IntegralType.interior_facet, subdomains) - - subdomains = self._subdomains.get("vertex") - if subdomains: - self._cpp_object.integrals.set_domains(cpp.fem.IntegralType.vertex, subdomains) + # Prepare dolfinx.cpp.fem.Form and hold it as a member + ffi = cffi.FFI() + self._cpp_object = cpp.fem.create_form(ffi.cast("uintptr_t", ufc_form), + function_spaces, coeffs, + [c._cpp_object for c in form.constants()], subdomains, mesh) diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp index fe2d08a182b..eaf323b2614 100644 --- a/python/dolfinx/wrappers/fem.cpp +++ b/python/dolfinx/wrappers/fem.cpp @@ -159,9 +159,17 @@ void fem(py::module& m) "create_form", [](const std::uintptr_t form, const std::vector< - std::shared_ptr>& spaces) { + std::shared_ptr>& spaces, + const std::vector>>& coefficients, + const std::vector>>& constants, + const std::map*>& subdomains, + const std::shared_ptr& mesh) { const ufc_form* p = reinterpret_cast(form); - return dolfinx::fem::create_form(*p, spaces); + return dolfinx::fem::create_form( + *p, spaces, coefficients, constants, subdomains, mesh); }, "Create Form from a pointer to ufc_form."); m.def( @@ -176,8 +184,6 @@ void fem(py::module& m) "build_dofmap", [](const MPICommWrapper comm, const dolfinx::mesh::Topology& topology, const dolfinx::fem::ElementDofLayout& element_dof_layout) { - // See https://github.com/pybind/pybind11/issues/1138 on why we need - // to convert from a std::unique_ptr to a std::shard_ptr auto [map, dofmap] = dolfinx::fem::DofMapBuilder::build( comm.get(), topology, element_dof_layout); return std::pair(map, std::move(dofmap)); @@ -194,7 +200,7 @@ void fem(py::module& m) .def(py::init([](const std::uintptr_t ufc_element) { const ufc_finite_element* p = reinterpret_cast(ufc_element); - return std::make_unique(*p); + return dolfinx::fem::FiniteElement(*p); })) .def("num_sub_elements", &dolfinx::fem::FiniteElement::num_sub_elements) .def("dof_reference_coordinates", @@ -353,83 +359,66 @@ void fem(py::module& m) // }, // py::return_value_policy::take_ownership); - // dolfinx::fem::FormIntegrals - py::class_, - std::shared_ptr>> - formintegrals(m, "FormIntegrals", - "Holder for integral kernels and domains"); - formintegrals - .def("integral_ids", - &dolfinx::fem::FormIntegrals::integral_ids) - .def_property_readonly( - "needs_permutation_data", - &dolfinx::fem::FormIntegrals::needs_permutation_data) - .def( - "integral_domains", - [](dolfinx::fem::FormIntegrals& self, - dolfinx::fem::IntegralType type, int i) { - const std::vector& domains - = self.integral_domains(type, i); - return py::array_t(domains.size(), domains.data(), - py::none()); - }, - py::return_value_policy::reference_internal, - "Return active domains for given integral") - .def("set_domains", - &dolfinx::fem::FormIntegrals::set_domains); - py::enum_(m, "IntegralType") .value("cell", dolfinx::fem::IntegralType::cell) .value("exterior_facet", dolfinx::fem::IntegralType::exterior_facet) - .value("interior_facet", dolfinx::fem::IntegralType::interior_facet); + .value("interior_facet", dolfinx::fem::IntegralType::interior_facet) + .value("vertex", dolfinx::fem::IntegralType::vertex); // dolfinx::fem::Form py::class_, std::shared_ptr>>( m, "Form", "Variational form object") - .def(py::init< - std::vector>, - bool>()) - .def_property_readonly("integrals", - &dolfinx::fem::Form::integrals) - .def_property_readonly( - "coefficients", - py::overload_cast<>(&dolfinx::fem::Form::coefficients, - py::const_)) - .def( - "num_coefficients", - [](const dolfinx::fem::Form& self) { - return self.coefficients().size(); - }, - "Return number of coefficients in form") - .def("original_coefficient_position", - [](dolfinx::fem::Form& self, int i) { - return self.coefficients().original_position(i); - }) - .def("set_coefficient", - [](dolfinx::fem::Form& self, std::size_t i, - std::shared_ptr> - f) { self.coefficients().set(i, f); }) - .def("set_constants", - [](dolfinx::fem::Form& self, - const std::vector>>& constants) { - auto& c = self.constants(); - if (constants.size() != c.size()) - throw std::runtime_error("Incorrect number of constants."); - for (std::size_t i = 0; i < constants.size(); ++i) - c[i] = std::pair("", constants[i]); - }) - .def("set_mesh", &dolfinx::fem::Form::set_mesh) - .def("set_tabulate_tensor", - [](dolfinx::fem::Form& self, - dolfinx::fem::IntegralType type, int i, py::object addr) { - auto tabulate_tensor_ptr = (void (*)( - PetscScalar*, const PetscScalar*, const PetscScalar*, - const double*, const int*, const std::uint8_t*, - const std::uint32_t))addr.cast(); - self.set_tabulate_tensor(type, i, tabulate_tensor_ptr); - }) + .def(py::init([](const std::vector>& spaces, + const std::map< + dolfinx::fem::IntegralType, + std::pair>, + const dolfinx::mesh::MeshTags*>>& + integrals, + const std::vector>>& + coefficients, + const std::vector>>& + constants, + bool needs_permutation_data, + const std::shared_ptr& mesh) { + using kern = std::function; + std::map>, + const dolfinx::mesh::MeshTags*>> + _integrals; + + // Loop over kernel for each entity type + for (auto& kernel_type : integrals) + { + // Set subdomain markers + _integrals[kernel_type.first].second = nullptr; + + // Loop over each domain kernel + for (auto& kernel : kernel_type.second.first) + { + auto tabulate_tensor_ptr = (void (*)( + PetscScalar*, const PetscScalar*, const PetscScalar*, + const double*, const int*, const std::uint8_t*, + const std::uint32_t))kernel.second.cast(); + _integrals[kernel_type.first].first.push_back( + {kernel.first, tabulate_tensor_ptr}); + } + } + return dolfinx::fem::Form( + spaces, _integrals, coefficients, constants, + needs_permutation_data, mesh); + }), + py::arg("spaces"), py::arg("integrals"), py::arg("coefficients"), + py::arg("constants"), py::arg("need_permutation_data"), + py::arg("mesh") = py::none()) + .def_property_readonly("coefficients", + &dolfinx::fem::Form::coefficients) .def_property_readonly("rank", &dolfinx::fem::Form::rank) .def_property_readonly("mesh", &dolfinx::fem::Form::mesh) .def_property_readonly("function_spaces", diff --git a/python/dolfinx/wrappers/mesh.cpp b/python/dolfinx/wrappers/mesh.cpp index 0ba3a02288d..a1ea62a0844 100644 --- a/python/dolfinx/wrappers/mesh.cpp +++ b/python/dolfinx/wrappers/mesh.cpp @@ -44,8 +44,8 @@ void declare_meshtags(py::module& m, std::string type) std::vector indices_vec(indices.data(), indices.data() + indices.size()); std::vector values_vec(values.data(), values.data() + values.size()); - return std::make_unique>( - mesh, dim, std::move(indices_vec), std::move(values_vec)); + return dolfinx::mesh::MeshTags(mesh, dim, std::move(indices_vec), + std::move(values_vec)); })) .def_readwrite("name", &dolfinx::mesh::MeshTags::name) .def_property_readonly("dim", &dolfinx::mesh::MeshTags::dim) @@ -173,7 +173,7 @@ void mesh(py::module& m) m, "Topology", "Topology object") .def(py::init([](const MPICommWrapper comm, const dolfinx::mesh::CellType cell_type) { - return std::make_unique(comm.get(), cell_type); + return dolfinx::mesh::Topology(comm.get(), cell_type); })) .def("set_connectivity", &dolfinx::mesh::Topology::set_connectivity) .def("set_index_map", &dolfinx::mesh::Topology::set_index_map) @@ -209,8 +209,7 @@ void mesh(py::module& m) .def(py::init([](const MPICommWrapper comm, const dolfinx::mesh::Topology& topology, dolfinx::mesh::Geometry& geometry) { - return std::make_unique(comm.get(), topology, - geometry); + return dolfinx::mesh::Mesh(comm.get(), topology, geometry); })) .def_property_readonly( "geometry", py::overload_cast<>(&dolfinx::mesh::Mesh::geometry), diff --git a/python/test/unit/fem/test_assemble_cppimport.py b/python/test/unit/fem/test_assemble_cppimport.py index 8182cb93aea..b4a7b0e125a 100644 --- a/python/test/unit/fem/test_assemble_cppimport.py +++ b/python/test/unit/fem/test_assemble_cppimport.py @@ -65,8 +65,8 @@ def compile_eigen_csr_assembler_module(): dolfinx::fem::assemble_matrix(mat_add, a, bcs); -auto map0 = a.function_space(0)->dofmap()->index_map; -auto map1 = a.function_space(1)->dofmap()->index_map; +auto map0 = a.function_spaces().at(0)->dofmap()->index_map; +auto map1 = a.function_spaces().at(1)->dofmap()->index_map; Eigen::SparseMatrix mat( map0->block_size() * (map0->size_local() + map0->num_ghosts()), map1->block_size() * (map1->size_local() + map1->num_ghosts())); diff --git a/python/test/unit/fem/test_custom_assembler.py b/python/test/unit/fem/test_custom_assembler.py index 335ae8b0cf0..47129b0ab29 100644 --- a/python/test/unit/fem/test_custom_assembler.py +++ b/python/test/unit/fem/test_custom_assembler.py @@ -127,7 +127,7 @@ def get_matsetvalues_api(): const PetscScalar* y, InsertMode addv); """) ffibuilder.set_source(module_name, """ - # include "petscmat.h" + #include "petscmat.h" """, libraries=['petsc'], include_dirs=[os.path.join(petsc_dir, petsc_arch, 'include'), diff --git a/python/test/unit/fem/test_custom_jit_kernels.py b/python/test/unit/fem/test_custom_jit_kernels.py index 1865dfb78e2..e7f74549445 100644 --- a/python/test/unit/fem/test_custom_jit_kernels.py +++ b/python/test/unit/fem/test_custom_jit_kernels.py @@ -74,13 +74,13 @@ def test_numba_assembly(): mesh = UnitSquareMesh(MPI.COMM_WORLD, 13, 13) V = FunctionSpace(mesh, ("Lagrange", 1)) - a = cpp.fem.Form([V._cpp_object, V._cpp_object], False) - a.set_tabulate_tensor(IntegralType.cell, -1, tabulate_tensor_A.address) - a.set_tabulate_tensor(IntegralType.cell, 12, tabulate_tensor_A.address) - a.set_tabulate_tensor(IntegralType.cell, 2, tabulate_tensor_A.address) + integrals = {IntegralType.cell: ([(-1, tabulate_tensor_A.address), + (12, tabulate_tensor_A.address), + (2, tabulate_tensor_A.address)], None)} + a = cpp.fem.Form([V._cpp_object, V._cpp_object], integrals, [], [], False) - L = cpp.fem.Form([V._cpp_object], False) - L.set_tabulate_tensor(IntegralType.cell, -1, tabulate_tensor_b.address) + integrals = {IntegralType.cell: ([(-1, tabulate_tensor_b.address)], None)} + L = cpp.fem.Form([V._cpp_object], integrals, [], [], False) A = dolfinx.fem.assemble_matrix(a) A.assemble() @@ -102,9 +102,8 @@ def test_coefficient(): vals = Function(DG0) vals.vector.set(2.0) - L = cpp.fem.Form([V._cpp_object], False) - L.set_tabulate_tensor(IntegralType.cell, -1, tabulate_tensor_b_coeff.address) - L.set_coefficient(0, vals._cpp_object) + integrals = {IntegralType.cell: ([(-1, tabulate_tensor_b_coeff.address)], None)} + L = cpp.fem.Form([V._cpp_object], integrals, [vals._cpp_object], [], False) b = dolfinx.fem.assemble_vector(L) b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) @@ -219,13 +218,13 @@ def test_cffi_assembly(): mesh.mpi_comm().Barrier() from _cffi_kernelA import ffi, lib - a = cpp.fem.Form([V._cpp_object, V._cpp_object], False) ptrA = ffi.cast("intptr_t", ffi.addressof(lib, "tabulate_tensor_poissonA")) - a.set_tabulate_tensor(IntegralType.cell, -1, ptrA) + integrals = {IntegralType.cell: ([(-1, ptrA)], None)} + a = cpp.fem.Form([V._cpp_object, V._cpp_object], integrals, [], [], False) - L = cpp.fem.Form([V._cpp_object], False) ptrL = ffi.cast("intptr_t", ffi.addressof(lib, "tabulate_tensor_poissonL")) - L.set_tabulate_tensor(IntegralType.cell, -1, ptrL) + integrals = {IntegralType.cell: ([(-1, ptrL)], None)} + L = cpp.fem.Form([V._cpp_object], integrals, [], [], False) A = dolfinx.fem.assemble_matrix(a) A.assemble() diff --git a/python/test/unit/fem/test_element_integrals.py b/python/test/unit/fem/test_element_integrals.py index 7d12b1812de..25c7cd02bca 100644 --- a/python/test/unit/fem/test_element_integrals.py +++ b/python/test/unit/fem/test_element_integrals.py @@ -128,8 +128,8 @@ def two_unit_cells(cell_type, agree=False, random_order=True, return_order=False return mesh -@ skip_in_parallel -@ parametrize_cell_types +@skip_in_parallel +@parametrize_cell_types def test_facet_integral(cell_type): """Test that the integral of a function over a facet is correct""" for count in range(5): @@ -168,8 +168,8 @@ def test_facet_integral(cell_type): assert np.isclose(result, out[0]) -@ skip_in_parallel -@ parametrize_cell_types +@skip_in_parallel +@parametrize_cell_types def test_facet_normals(cell_type): """Test that FacetNormal is outward facing""" for count in range(5): @@ -229,9 +229,9 @@ def test_facet_normals(cell_type): assert ones == 1 -@ skip_in_parallel -@ pytest.mark.parametrize('space_type', ["CG", "DG"]) -@ parametrize_cell_types +@skip_in_parallel +@pytest.mark.parametrize('space_type', ["CG", "DG"]) +@parametrize_cell_types def test_plus_minus(cell_type, space_type): """Test that ('+') and ('-') give the same value for continuous functions""" results = [] @@ -249,9 +249,9 @@ def test_plus_minus(cell_type, space_type): assert np.isclose(i, j) -@ skip_in_parallel -@ pytest.mark.parametrize('pm', ["+", "-"]) -@ parametrize_cell_types +@skip_in_parallel +@pytest.mark.parametrize('pm', ["+", "-"]) +@parametrize_cell_types def test_plus_minus_simple_vector(cell_type, pm): """Test that ('+') and ('-') match up with the correct DOFs for DG functions""" results = [] @@ -299,10 +299,10 @@ def test_plus_minus_simple_vector(cell_type, pm): assert np.isclose(results[0][dof0], result[dof1]) -@ skip_in_parallel -@ pytest.mark.parametrize('pm1', ["+", "-"]) -@ pytest.mark.parametrize('pm2', ["+", "-"]) -@ parametrize_cell_types +@skip_in_parallel +@pytest.mark.parametrize('pm1', ["+", "-"]) +@pytest.mark.parametrize('pm2', ["+", "-"]) +@parametrize_cell_types def test_plus_minus_vector(cell_type, pm1, pm2): """Test that ('+') and ('-') match up with the correct DOFs for DG functions""" results = [] @@ -353,10 +353,10 @@ def test_plus_minus_vector(cell_type, pm1, pm2): assert np.isclose(results[0][dof0], result[dof1]) -@ skip_in_parallel -@ pytest.mark.parametrize('pm1', ["+", "-"]) -@ pytest.mark.parametrize('pm2', ["+", "-"]) -@ parametrize_cell_types +@skip_in_parallel +@pytest.mark.parametrize('pm1', ["+", "-"]) +@pytest.mark.parametrize('pm2', ["+", "-"]) +@parametrize_cell_types def test_plus_minus_matrix(cell_type, pm1, pm2): """Test that ('+') and ('-') match up with the correct DOFs for DG functions""" results = [] @@ -409,9 +409,9 @@ def test_plus_minus_matrix(cell_type, pm1, pm2): assert np.isclose(results[0][a, c], result[b, d]) -@ skip_in_parallel -@ pytest.mark.parametrize('order', [1, 2]) -@ pytest.mark.parametrize('space_type', ["N1curl", "N2curl"]) +@skip_in_parallel +@pytest.mark.parametrize('order', [1, 2]) +@pytest.mark.parametrize('space_type', ["N1curl", "N2curl"]) def test_curl(space_type, order): """Test that curl is consistent for different cell permutations of a tetrahedron."""