How do I input gridwise data?

General
1

Hello All,

I am using electric and magnetic fields as initialized FEM data in my simulation, and I am shifting from using analytical definitions to applying numerically calculated fields (i.e. a seperate function that outputs a list of E and B fields for each grid point).

How do I then input the fields (now numpy arrays of E1, E2, E3, B1, B2, B3 over each grid node) into getFEM as initialized FEM data? Is this possible?

Sincerely,
Eric Comstock

2

Maybe you can find the answer by studying my maze solver code:

image
image512×288 90.4 KB

The maze is read from a bitmap file as grayscale values in a fixed grid, and interpolated to a finite element mesh.

import getfem as gf
import numpy as np
import PIL.Image,PIL.ImageOps

img = PIL.Image.open("maze.png")
data = np.array(img)
data[:,:,:][data[:,:,3]<255] = [255,255,255,255]
img = PIL.ImageOps.grayscale(PIL.Image.fromarray(data))
data = np.array(img)

nonemptyrows = np.where(np.any(data[:,:]!=255,axis=1))[0]
nonemptycols = np.where(np.any(data[:,:]!=255,axis=0))[0]
data = data[nonemptyrows[0]:nonemptyrows[-1],nonemptycols[0]:nonemptycols[-1]] // 255

NY,NX = data.shape    # number of elements in horizontal and vertical direction

B_RG = 99 # bottom mesh region
T_RG = 98 # top mesh region
m = gf.Mesh("import", "structured",
            f"GT='GT_QK(2,1)';ORG=[-0.5,-0.5];SIZES=[{NX},{NY}];NSUBDIV=[{NX//4},{NY//4}]")
m.set_region(B_RG, m.outer_faces_with_direction([0,-1], np.pi/10))
m.set_region(T_RG, m.outer_faces_with_direction([0,1], np.pi/10))

mfu = gf.MeshFem(m,2)    # vector fem with 2 components
mfu.set_classical_fem(2) # 2nd order Lagrange fem (9-node quads or 6-node trias)

mfp = gf.MeshFem(m,1)    # scalar fem
mfp.set_classical_fem(1) # 1st order Lagrange fem (4-node quads or 3-node trias)

mfrho = gf.MeshFem(m,1)    # scalar fem
mfrho.set_classical_fem(0) # elementwise constant fem (discontinuous)

mim = gf.MeshIm(m, 5)

md = gf.Model("real")
md.add_fem_variable("u", mfu)
md.add_fem_variable("p", mfp)

md.add_fem_data("rho", mfrho)
X = md.interpolation("X(1)", mfrho).round().astype(int)
Y = md.interpolation("X(2)", mfrho).round().astype(int)
md.set_variable("rho", data[NY-Y-1,X])

md.add_linear_term(mim, "(1-0.99999999*rho)*1e6*u.Test_u+Grad(u):Grad(Test_u)"
                        "+Grad(p).Test_u+Div(u)*Test_p") # stokes equations
md.add_linear_term(mim, "(p-100)*Test_p", T_RG)
md.add_linear_term(mim, "(p-0)*Test_p", B_RG)

md.solve("lsolver", "mumps")
mfrho.export_to_vtu("maze.vtu", mfrho, md.variable("rho"), "rho",
                                mfu, md.variable("u"), "Velocity",
                                mfp, md.variable("p"), "Pressure")
3

here you can find another example

https://gitlab.gbar.dtu.dk/-/snippets/44

repeating the snippet her for your convenience:

import getfem as gf
import cv2

LX = ...
LY = ...
NX = ...

...

mfrho = ...

im = cv2.imread("./field.png", cv2.IMREAD_GRAYSCALE)
im = np.flip(cv2.GaussianBlur(im,(49,49),0),0) # mirror upside down and blur
im = im - np.min(im)
im = 1. - im/np.max(im)

NX0 = im.shape[1]
NY0 = im.shape[0]
m_aux = gf.Mesh("import", "structured",
                "GT='GT_QK(2,1)';ORG=[0,0];SIZES=[%e,%e];NSUBDIV=[%i,%i]"
                % (LX, LY, NX0, NY0))
mf_aux = gf.MeshFem(m_aux, 1)
mf_aux.set_classical_fem(0)
mim_aux = gf.MeshIm(m_aux, 3)

md_aux = gf.Model("real")
md_aux.add_fem_data("rho0", mf_aux)
md_aux.set_variable("rho0", np.reshape(im, NY0))
md_aux.add_fem_variable("rho", mfrho)
md_aux.add_interpolate_transformation_from_expression("mesh1to2", m, m_aux, "X")
md_aux.add_linear_term(mim, "(rho_tilde-Interpolate(rho0,mesh1to2))*Test_rho_tilde")
md_aux.solve()
RHO = md_aux.variable("rho")
4

Hello Kostas,

Thank you so much! I will see if I can figure it out.

Sincerely,

Eric Comstock