Calculate cell data

[martinbaun]

Hi,

I`m running a simulation of an Iron Core with a current density driven coil. The mesh is shown below. The attached code leads to correct solution vector "A". For postprocessing I want to calculate the magnetic flux density. In ex35 there is an example how to do that as nodal values. This leads to a blurry image and the image shows a magnetic field outside the iron core which shouldn't be. So I would like to calulate the flux density as a cell value. For linear elements the flux density is a constant in the cell, which is shown in the last picture below which I created with an other software. Do you have any hints?

Best regards
Martin

import numpy as np
mu_0 = np.pi*4e-7
import matplotlib.pyplot as plt
plt.ion()
import skfem
from skfem.helpers import dot, grad


m = skfem.Mesh.load('mesh.msh')
e = skfem.ElementTriP1()
basis = skfem.Basis(m, e)

def PolyArea(x,y):
    return 0.5*np.abs(np.dot(x,np.roll(y,1))-np.dot(y,np.roll(x,1)))
    
area = []
for k in range(m.nelements):
    nodes = m.t[:,k]
    x = m.p[0,:][nodes]
    y = m.p[1,:][nodes]
    area.append(PolyArea(x, y))


perm = np.ones(m.nelements)
perm[m.subdomains['Anke']] = 1000
perm[m.subdomains['E-Ke']] = 1000
perm *= mu_0
perm = np.array([perm, perm, perm]).T

# current density
currdens_element = np.zeros(m.nelements)
currdens_element[m.subdomains['Winding_1_+']] = 10*100/(8e-3*13e-3)
currdens_element[m.subdomains['Winding_1_-']] = -10*100/(8e-3*13e-3)

J = np.zeros(m.nvertices)
for k in range(len(currdens_element)):
    nodes = m.t[:,k]
    J[nodes] += currdens_element[k]/len(nodes)*area[k]


@skfem.BilinearForm
def laplace(u, v, _):
    return dot(grad(u), grad(v))/perm


# this method could also be imported from skfem.models.unit_load
@skfem.LinearForm
def rhs(v, _):
    return 1.0 * v


S = skfem.asm(laplace, basis)
b = J
S, b = skfem.enforce(S, b, D=m.boundary_nodes())
A = skfem.solve(S, b)

B_x = skfem.projection(A, basis, basis, 1)
B_y = -skfem.projection(A, basis, basis, 0)

from skfem.visuals.matplotlib import plot
plot(basis, np.sqrt(B_x**2 + B_y**2), colorbar=True)
plt.axis('equal')

[kinnala]

You could try creating another basis object with piecewise-constant element and project to that:

basis0 = basis.with_element(skfem.ElementTriP0())
B_x = skfem.projection(A, basis0, basis, 1)
B_y = -skfem.projection(A, basis0, basis, 0)

from skfem.visuals.matplotlib import plot
plot(basis0, np.sqrt(B_x**2 + B_y**2), colorbar=True)
plt.axis('equal')

Let me know if it works or if there is another issue.

[martinbaun]

Wow, that was a fast response. Works perfect!