### Multimesh with different material properties for each part

153
views
0
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
{
public:
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));
}
public:
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 'ft11_magnetostatics.py'--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
shims = [generate_mesh(Circle(Point(rPos*cos(v), rPos*sin(v)), shimRad), 10) for v in pos]

mm = MultiMesh()
for item in shims:
mm.build()

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:
W.build()

h = MultiMeshFunction(W)



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

h.assign_part(0,Constant(0))

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.

Cheers!

Community: FEniCS Project

2
9 weeks ago by
Hi,

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