Add basic finite element matrices construction in python

feelpp/pyfeelpp/feelpp/CMakeLists.txt

@@ -12,12 +12,13 @@ feelpp_add_pymodule(vf SRCS vf.cpp DESTINATION feelpp )
 feelpp_add_pymodule(measure SRCS measure.cpp DESTINATION feelpp )
 feelpp_add_pymodule(integrate SRCS integrate.cpp DESTINATION feelpp )
 feelpp_add_pymodule(interpolation SRCS interpolation.cpp DESTINATION feelpp )
+feelpp_add_pymodule(operators SRCS operators.cpp DESTINATION feelpp )
 feelpp_add_pymodule(models SRCS models.cpp DESTINATION feelpp )
 feelpp_add_pymodule(quality SRCS quality.cpp DESTINATION feelpp )
 
-foreach(MODULE IN ITEMS core feelpp alg mesh discr exporter ls ts vf measure models integrate interpolation quality)
+foreach(MODULE IN ITEMS core feelpp alg mesh discr exporter ls ts vf measure models integrate interpolation operators quality)
   add_dependencies(pyfeelpp _${MODULE})
 endforeach()
 
-set(PYFILES __init__.py integrate.py measure.py quality.py interpolation.py )
+set(PYFILES __init__.py integrate.py measure.py quality.py interpolation.py operators.py )
 install(FILES ${PYFILES} DESTINATION ${FEELPP_PYTHON_MODULE_PATH}/feelpp)

feelpp/pyfeelpp/feelpp/operators.cpp

@@ -0,0 +1,110 @@
+//! -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t  -*- vim:fenc=utf-8:ft=cpp:et:sw=4:ts=4:sts=4
+//!
+//! This file is part of the Feel++ library
+//!
+//! This library is free software; you can redistribute it and/or
+//! modify it under the terms of the GNU Lesser General Public
+//! License as published by the Free Software Foundation; either
+//! version 2.1 of the License, or (at your option) any later version.
+//!
+//! This library is distributed in the hope that it will be useful,
+//! but WITHOUT ANY WARRANTY; without even the implied warranty of
+//! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+//! Lesser General Public License for more details.
+//!
+//! You should have received a copy of the GNU Lesser General Public
+//! License along with this library; if not, write to the Free Software
+//! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+//!
+//! @file
+//! @author Christophe Prud'homme <christophe.prudhomme@feelpp.org>
+//! @date 23 Jul 2021
+//! @copyright 2017 Feel++ Consortium
+//!
+#include <pybind11/pybind11.h>
+#include <fmt/core.h>
+#include <feel/feelcore/environment.hpp>
+#include <feel/feeldiscr/mesh.hpp>
+#include <feel/feeldiscr/operatorinterpolation.hpp>
+#include <feel/feelfilters/loadmesh.hpp>
+#include <feel/feelmesh/filters.hpp>
+#include <feel/feelvf/vf.hpp>
+#include <mpi4py/mpi4py.h>
+#include <pybind11/eigen.h>
+#include <boost/algorithm/string.hpp>
+namespace py = pybind11;
+
+template <typename SpaceT>
+void defOperator( py::module& m )
+{
+    using namespace Feel;
+    using space_t = SpaceT;
+    using space_ptr_t = std::shared_ptr<space_t>;
+    using mesh_t = typename SpaceT::mesh_type;
+    using mesh_ptr_t = std::shared_ptr<mesh_t>;
+    constexpr int Dim = mesh_t::nDim;
+    constexpr int Order = space_t::basis_type::nOrder;
+    constexpr int RealDim = mesh_t::nRealDim;
+    using size_type = uint32_type;
+
+
+    std::string suffix;
+    if ( space_t::is_continuous && space_t::is_scalar )
+        suffix = std::string("Pch");
+    if ( space_t::is_continuous && space_t::is_vectorial )
+        suffix = std::string("Pchv");
+    if ( !space_t::is_continuous && space_t::is_scalar )
+        suffix = std::string("Pdh");
+    if ( !space_t::is_continuous && space_t::is_vectorial )
+        suffix = std::string("Pdhv");
+    std::string pyclass_name = fmt::format( "Mass_{}_{}D_P{}", suffix, Dim, Order );
+    VLOG(2) << fmt::format("[pyfeelpp] class name: {}", pyclass_name ) << std::endl;
+    using iterator_range_t = elements_reference_wrapper_t<mesh_t>;
+
+    m.def(
+        "mass", []( space_ptr_t const& domain, space_ptr_t const& image, iterator_range_t const& r, std::string const& c = "1" )
+        { 
+            auto a=form2(_test=image,_trial=domain);
+            auto u=domain->element();
+            auto v=image->element();
+
+            a=integrate(_range=r,_expr=expr(c)*trace(trans(idt(u))*id(v)));
+            return a.matrixPtr(); },
+        py::arg( "trial" ), py::arg( "test" ), py::arg( "range" ), py::arg( "coeff" )=std::string{"1"}, "create a mass matrix" );
+    m.def(
+        "stiffness", []( space_ptr_t const& domain, space_ptr_t const& image, iterator_range_t const& r, std::string const& c = "1" )
+        { 
+            auto a=form2(_test=image,_trial=domain);
+            auto u=domain->element();
+            auto v=image->element();
+            if ( boost::algorithm::contains(c, "{") && boost::algorithm::contains(c, "}") )
+                a = integrate( _range = r, _expr = trace( expr<RealDim, RealDim>( c ) * trans( gradt( u ) ) * grad( v ) ) );
+            else
+                a=integrate(_range=r,_expr=trace(expr(c)*trans(gradt(u))*grad(v)));
+            return a.matrixPtr(); },
+        py::arg( "trial" ), py::arg( "test" ), py::arg( "range" ), py::arg( "coeff" )=std::string{"1"},"create a stiffness matrix" );
+}
+
+PYBIND11_MODULE( _operators, m )
+{
+    using namespace Feel;
+    using namespace hana::literals;
+    if ( import_mpi4py() < 0 ) return;
+
+    auto ordert = hana::make_tuple( 1_c, 2_c );
+    hana::for_each( ordert, [&m](auto const& o ){
+        constexpr int _order = std::decay_t<decltype(o)>::value;
+        // 1D
+        //std::cout << fmt::format("-- Pch 1D P{}", _order ) << std::endl;
+        defOperator<Pch_type<Mesh<Simplex<1>>, _order>>( m );
+        // 2D
+        //std::cout << fmt::format("-- Pch 2D P{}", _order ) << std::endl;
+        defOperator<Pch_type<Mesh<Simplex<2>>, _order>>( m );
+        //std::cout << fmt::format("-- Pchv 2D P{}", _order ) << std::endl;
+        defOperator<Pchv_type<Mesh<Simplex<2>>, _order>>( m );
+        // 3D
+        //std::cout << fmt::format("-- Pch 3D P{}", _order ) << std::endl;
+        defOperator<Pch_type<Mesh<Simplex<3>>, _order>>( m );
+    });
+
+}

feelpp/pyfeelpp/feelpp/operators.py

@@ -0,0 +1,3 @@
+from ._operators import *
+
+

feelpp/pyfeelpp/tests/test_operators.py

@@ -0,0 +1,31 @@
+import math
+import sys
+import feelpp
+from feelpp.operators import *
+import pytest
+
+
+def run(m, geo):
+    mesh_name, e_meas, e_s_1, e_s_2, e_s_bdy=geo
+    mesh= feelpp.load(m, mesh_name, 0.1)
+    Xh = feelpp.functionSpace(mesh=mesh)
+    v=Xh.element()
+    v.on(range=feelpp.elements(mesh), expr=feelpp.expr("1"))
+
+    M=mass(test=Xh,trial=Xh,range=feelpp.elements(mesh))
+    assert(abs(M.energy(v,v)-e_meas)<1e-10)
+    S=stiffness(test=Xh,trial=Xh,range=feelpp.elements(mesh))
+    assert(abs(S.energy(v,v)-0)<1e-10)
+    v.on(range=feelpp.elements(mesh), expr=feelpp.expr("x+y:x:y"))
+    assert(abs(S.energy(v,v)-2*e_meas)<1e-10)
+
+
+def test_operators(init_feelpp):
+    feelpp.Environment.changeRepository(
+        directory="pyfeelpp-tests/measure")
+    geo = {
+        '2': feelpp.create_rectangle(),
+        '3': feelpp.create_box(),
+    }
+    run(feelpp.mesh(dim=2), geo['2'])
+    run(feelpp.mesh(dim=3, realdim=3), geo['3'])

pyfeelpp-all/CMakeLists.txt

@@ -1,4 +1,4 @@
-cmake_minimum_required(VERSION 3.3)
+cmake_minimum_required(VERSION 3.21)
 project(pyfeelpp-all)
 
 if ( (CMAKE_SOURCE_DIR STREQUAL CMAKE_CURRENT_SOURCE_DIR) OR (FEELPP_COMPONENT STREQUAL "pyfeelpp-all" ))