### How to access solution values in c++ coded expression?

344
views
0
7 months ago by
Hi,
I'm trying to solve a nonlinear heat equation, with thermal conductivity depending on the material temperature. I have tabulated data for the conductivity and want to interpolate it in a C++ coded Expression, so it works fast. However, I do not know how to access the solution value from within that C++ code snippet. Consider the following example code:
from fenics import *

# Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
u_D = Constant(0.0)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, u_D, boundary)

# Define variational problem
u = Function(V)
du = TrialFunction(V)
v = TestFunction(V)
f = Constant(6.0)

interpolatedConductivity = Expression(
'''
// example of linear interpolation of data points

class interpolatedConductivity : public Expression
{
public:

const double kdata[7] = {1,1,2,3,4,5,5}; // an array of data points corresponding to temperatures at {0, 5, 10, .. 30}
const double minT = 0.0;
const double maxT = 30.0;
const double dT = 5.0;

interpolatedConductivity() : Expression()
{
}

// This eval function is passed some arguments, how can I get the solution value as an input argument? So x would be the temperature in an element.
void eval(double& value, const double& x) const
{
if(x > maxT)            // do some check whether input is in table
value = kdata[6];
else if(x < minT)
value = kdata[0];
else
{
int lowerIndex = floor(x - minT);
int upperIndex = ceil(x - minT);
double relT = x - minT - lowerIndex*dT; // relative position between the two data points
value = (1 - relT)*kdata[lowerIndex] + relT*kdata[upperIndex]; // linear interpolation
}
}
};
''', degree = 1)

F = interpolatedConductivity(u)*dot(grad(u), grad(v))*dx - f*v*dx # does not work at all

J = derivative(F, u, du)
problem = NonlinearVariationalProblem(F, u, bc, J=J)
solver = NonlinearVariationalSolver(problem)

# Compute solution
solver.solve()
​
I know Expressions are usually passed a coordinate vector as default. What could I change n the declaration of "void eval..." to let it access the current solution vector?
Or is this a completely wrong approach?

Thanks in advance for any help or hint!
Community: FEniCS Project

1
7 months ago by
Here's a very ugly solution using Python, but sticking to UFL, so it gets compiled into C++:

from fenics import *

# Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
u_D = Constant(0.0)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, u_D, boundary)

# Define variational problem
u = Function(V)
du = TrialFunction(V)
v = TestFunction(V)
f = Constant(6.0)

k0 = 1.0
k1 = 1.0
k2 = 2.0
k3 = 3.0
k4 = 4.0
k5 = 5.0
k6 = 5.0

t0 = 0.0
t1 = 5.0
t2 = 10.0
t3 = 15.0
t4 = 20.0
t5 = 25.0
t6 = 30.0

def c(a,b,c):
return conditional(a,b,c)

def interp(t,thi,tlo,khi,klo):
frac = (t-tlo)/(thi-tlo)
return frac*khi + (1.0-frac)*klo

kappa = c(gt(u,t0),
c(gt(u,t1),
c(gt(u,t2),
c(gt(u,t3),
c(gt(u,t4),
c(gt(u,t5),
c(gt(u,t6),
k6,
interp(u,t5,t6,k5,k6)),
interp(u,t4,t5,k4,k5)),
interp(u,t3,t4,k3,k4)),
interp(u,t2,t3,k2,k3)),
interp(u,t1,t2,k1,k2)),
interp(u,t0,t1,k0,k1)),
k0)

J = derivative(F, u, du)
problem = NonlinearVariationalProblem(F, u, bc, J=J)
solver = NonlinearVariationalSolver(problem)

# Compute solution
solver.solve()
​

You should double-check my logic, but it does run, and hopefully the idea is clear.  I doubt this sort of conditional() abuse is considered best practice, though, especially if your real data has many more temperature/conductivity pairs...
Yes, I also had that idea, but there are indeed hundreds of measurement points ( and I would do that for multiple coefficients ), so I assume it would be annoyingly slow in the end. Thank you though for contributing!
written 7 months ago by julianL
Okay, this piqued my interest, so I dug up a way to call custom C++ code that involves DOLFIN objects.  My starting point was the code in hippylib

https://github.com/hippylib/hippylib

for assembling pointwise observations.  (I am not involved in that project; I just happened to remember that it had a function like this.)  I put together a minimal example, involving three files: a Python script, "p.py"

from dolfin import *
import os

mesh = UnitSquareMesh(10,10)
VE = FiniteElement("Lagrange",mesh.ufl_cell(),1)
V = FunctionSpace(mesh,VE)

f = Function(V)

# Giving credit: I've paraphrased Umberto Villa's code
#
# https://github.com/hippylib/hippylib/blob/master/hippylib/pointwiseObservation.py
#
# (and associated C++ files)
#
# to figure out how to compile custom C functions that use DOLFIN objects.
# (I'd been meaning to go through it for some time to figure out how to assemble
# things at arbitrary points, to implement immersed methods.)

abspath = os.path.dirname( os.path.abspath(__file__) )
source_directory = os.path.join(abspath,".")
cpp_sources = ["c.cpp",]

include_dirs = [".",source_directory]

# This doesn't work for me, although it's probably a quirk of my mixed-up
# environment variables.
include_dirs = [".", source_directory]
for ss in ['PROFILE_INSTALL_DIR', 'PETSC_DIR', 'SLEPC_DIR']:
if ss in os.environ.keys():
include_dirs.append(os.environ[ss]+'/include')

# I've just hard-coded my PETSc directory instead. This is specific to my
# Ubuntu installation.
include_dirs.append("/usr/lib/petscdir/3.7/include/")

cpp_module = compile_extension_module(
code = code, source_directory = source_directory,
sources = cpp_sources, include_dirs = include_dirs)

# The Function is, by default, zero, so this just prints zero; it's not
# very exciting, but the real difficulty here is getting everything to
# compile and run.
print(cpp_module.eval_f(f,Point(0.5,0.5)))

​​

#include <dolfin/function/Function.h>
#include <dolfin/geometry/Point.h>

namespace dolfin
{
double eval_f(Function &f, Point &x);
}
​

and a C++ source file, "c.cpp",

#include <dolfin/common/version.h>
#include <dolfin/common/Array.h>

#include "c.h"

using namespace dolfin;

double dolfin::eval_f(Function &f, Point &x){
Array<double> values(1);
Array<double> coords(2);
for(int i=0; i<2; i++)
coords[i] = x.coordinates()[i];
f.eval(values,coords);
return values[0];
} // end eval_f
​​

All of these go in one directory.  I had to hard-code a library path in "p.py", so read the comments and maybe adjust accordingly.  You may need to tinker a bit to get everything working on your platform.  (I have 2017.1 installed via the Ubuntu package.)  For getting an Expression, the most direct extension of this example would be to just create a sub-class of Expression in Python, then define its eval_cell() method to call your custom C++ function. You may also be able to directly create a subclass in your C++ module, but you'll probably need to fight your way through a lengthy series of cryptic errors (like I did to get this example working...).
written 7 months ago by David Kamensky
0
7 months ago by