How to create piecewise expression based on value of a spatially varying function?

variational formulation
1

Hello all,

Consider if I have a spatially varying function, A on some mesh, for example:

mesh = UnitSquareMesh(10,10)
V = FunctionSpace(mesh, 'Lagrange',1)
A = Function(V)

A in reality would be the output from a linear solve.
How would I implement a second function, B which is piecewise dependant on the values in A? i.e.:

B(A) =\begin{cases} B_0(A) & A < 1 \\ B_1(A) & 1 \leq A < 5 \\ B_2(A) & 5 \leq A < 10 \\ B_3(A) & A > 10 \end{cases}

Where B_0 \dots B_3 are functions in A, say, as an expression: B(A) = A^2 + 5 etc.

Both A and B would then require to appear in the variational formulation of my problem.

Any help would be greatly appreciated.

2

I’d solve it similar to subdomains as in Subdomains and boundary conditions : Use an Expression with the variable A as input (as an Expression) and then a simple if-loop.
Something like this shuld do the trick:

cpp_code =''' #include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
namespace py = pybind11;
#include <math.h> // to use pow()
#include <dolfin/function/Expression.h> // for expressions
#include <dolfin/mesh/MeshFunction.h> //  for MeshFunctions
#include <dolfin/function/Function.h> // for functions



class your_Function_Name : public dolfin::Expression
{
   public:
       your_Function_Name() : dolfin::Expression(){} 
       void eval(Eigen::Ref<Eigen::VectorXd> values, Eigen::Ref<const Eigen::VectorXd> x, const ufc::cell& cell) const override
       {

           exp1->eval(values, x);
           const double val1 = values[0];

          if ( val1 < 1){
                  values[0] = 10; // your stuff
          }else{
                 values[0] = pow(val1, 2)-3; // your stuff
          {

       }


       std::shared_ptr<dolfin::Expression> exp1;
};
PYBIND11_MODULE(SIGNATURE, m)
{
   py::class_<your_Function_Name, std::shared_ptr<your_Function_Name>, dolfin::Expression>(m, "your_Function_Name")
   .def(py::init<>())
   .def_readwrite("exp1", &your_Function_Name::exp1);
}
'''
B=CompiledExpression( compile_cpp_code(cpp_code).your_Function_Name(), exp1=A, degree=2)

Note that I haven’t tested it, but this should give you hopefully an idea :slight_smile:
Hope this helps :slight_smile:

3

Hello @Emanuel_Gebauer,

Thank you so much for your reply! - I should’ve specified, I’m unfamiliar with C++ and using Python. However, you pointed me in the right direction, this is what I came up with and it seems to work:

from fenics import *

# Function class
class piecewise_fcn(UserExpression):
    def set_values(self,A):
        self.A = A
    def eval(self, value, x):
        A_eval = A(x)
        if A_eval < 1:
            value[0] = 1
        elif A_eval >= 1 and A_eval < 5:
            value[0] = 2
        elif A_eval >= 5 and A < 10:
            value[0] = 3
        else:
            value[0] = 4
    def value_shape(self):
        return()

Please let me know if there is any way I could improve upon this, but i’ll leave this here as the solution for now,

Thanks again!