Multimesh with different material properties for each part

9 weeks ago by
I am trying to calculate the magnetic field from an array of permanent magnet elements by direct integration using assemble. I was able to use this technique for the case of a single magnetic source by defining an expression that depends on 'rx,ry' (point at which I want to calculate the field), and 'mx,my', magnetization components of mesh elements.

cpp_code = """
class MyCppExpression : public Expression
    MyCppExpression() : Expression(2), rx(0), ry(0), mx(0), my(0){ }
    void eval(Array<double>& values, const Array<double>& x) const
        const double dx = x[0] - rx; 
        const double dy = x[1] - ry;
        const double rmag = sqrt(dx*dx + dy*dy);
        const double mdotr = dx*mx + dy*my;
        const double mu0 = 4*pi*1e-7;
        values[0] = mu0/(4*pi) * ((4*dx*mdotr)/pow(rmag,4) - 2*mx/pow(rmag,2));
        values[1] = mu0/(4*pi) * ((4*dy*mdotr)/pow(rmag,4) - 2*my/pow(rmag,2));
    double rx, ry, mx, my;

Expanding this problem to multiple source with identical magnetization vectors is straightforward using MultiMesh. The difficulty I am having is assigning a different magnetization vector to each element of a given part. I tried following the example ''--using a meshfunction to assign material properties, but I am having difficulty understanding how to assign a fixed value to a part of a multimesh function.

nShims = 16
pos = [i*2*pi/nShims for i in range(nShims)]

rPos = 0.25
shimRad = 0.01
shims = [generate_mesh(Circle(Point(rPos*cos(v), rPos*sin(v)), shimRad), 10) for v in pos]

mm = MultiMesh()
for item in shims:

fp1 = Point(1,0)
nPts = 100
fp = [i*2*pi/nPts for i in range(nPts)]

# Define mm function space
W = MultiMeshFunctionSpace(mm,'CG',1)
for item in shims:

h = MultiMeshFunction(W)

When I try to assign a value to the a part multimeshfunction I get an error


TypeError: in method 'MultiMeshFunction_assign_part', argument 3 of type 'dolfin::Function const &'

The value I want to assign is an identifier I can use in a switch/case statement in my cpp_function.

Any advice on how to implement this is appreciated.


Community: FEniCS Project

1 Answer

9 weeks ago by

I don't know if it's the proper way to do it but try with
h.assign_part(0, interpolate(Constant(0), W._parts[0]))​
Thank you. This seems to work as expected.

I have two additional questions. Is it possible to assign a vector value to this function space?

How do I evaluate the function to see if it is behaving properly before integrating it into my cpp_code for a given ufc::cell or dolfin::cell?

written 9 weeks ago by mike  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »