Best practice for a combined integration and interpolation

General
1

Hello,

I would like to know if there are current best practices are for combined integration and interpolation on a solution.

For my example, I have a 6D mesh with a particle density on it, with three dimensions of space and three of momentum. I would like to integrate over the momentum (p) coordinates to get the density (integral density dp) and average velocity ((integral density*p/m dp) / (integral density dp)) at each point on the 3D coordinates of position to help with troubleshooting my code.

Is there a generally accepted method to do this? If not, would an integration and then an interpolation, or an interpolation and then an integration make the most sense?

Sincerely,
Eric Comstock

2

Hello Eric, I do not understand the question entirely because 6 dimensions are not very intuitive for us common people. But I guess what you mean.

When you have a 6D mesh, the volume is 6D and surfaces are 5D (regions defined by elm, number and face number). If you define an integration method on a 6D mesh in getfem, it will define integration points in the interior (6D) and on the faces of the elements (5D)

If you want to integrate on any entity below 5D, then you need to create a new mesh for that. You will need to create a mesh with 3D elements inside the 6D space and interpolate 6D fields to the Gauss points of the 3D elements, then you can do the integration on this second mesh.

The 3D equivalent is to define a line mesh in 3D space, then interpolate from some 3D volume mesh onto the line mesh, and then integrate along the line mesh.

It should be possible. To transfer the field from one mesh to the other you need to define the identity transformation like in the following example from my Elastohydrodynamic model:

md.add_interpolate_transformation_from_expression("surface", m2, m1, "X")
md.add_interpolate_transformation_from_expression("substrate", m1, m2, "X")
md.add_macro("H", "H0+sqr(X(1))/2-Interpolate(U,substrate)(2)")
1 Like
3

Hello Kostas,

If it makes it clearer, imagine I have a function f(x,y,z,u,v,w) in 6D.

Then, I want to get the integral of f du dv dw, which is itself still a function of x, y, and z. So therefore, over the 6D volume, three dimensions are integrated and the other three remain. I then want to interpolate this over a grid of points, or just directly extract the integration.

Does that make sense?

Sincerely,
Eric Comstock

4

yes, then I understood it correctly, study the following snippet. The 6D volume I will call V6, any hypersurfaces in dimensions 5,4,3,2… I will call S5, S4, S3, …

mV6 = gf.Mesh('regular', grid_x, grid_y, grid_z, grid_u, grid_v, grid_w)
mS3 = gf.Mesh('regular', grid_u,grid_v,grid_w)
T_3to6 = np.zeros((6,3))
T_3to6[3,0] = T_3to6[4,1] = T_3to6[5,2] = 1.
mS3.transform(T_3to6)

mfV6 = gf.MeshFem....

mimS3 = gf.MeshIm(mS3) 

md = gf.Model("real")
md.add_initialized_fem_data("somevariable", mfV6, SOMEVARIABLENODALDATA)
md.add_interpolate_transformation_from_expression("V6toS3", mS3, mV6, "X")

I_S3 = gf.asm_generic(mimS3, 0, "Interpolate(somevariable,V6toS3)", -1, md)

I_S3 is the integral you are looking for.

1 Like
5

Hello Kostas,

When I run the code with the md.add_initialized_fem_data("f", mf, density) line, it says that f is already defined

RuntimeError: (Getfem::InterfaceError) -- Error in getfem_models.cc, line 129 : 
Variable f already exists


logic_error exception caught

But when I run it without the line, it says that the variable f is not defined

RuntimeError: (Getfem::InterfaceError) -- Error in getfem_generic_assembly_compile_and_exec.cc, line 5659 void getfem::ga_compile_node(getfem::pga_tree_node, getfem::ga_workspace&, getfem::ga_instruction_set&, getfem::ga_instruction_set::region_mim_instructions&, const getfem::mesh&, bool, getfem::ga_if_hierarchy&): 
The finite element of variable f has to be defined on the same mesh than the integration method or interpolation used


logic_error exception caught

Right now I am trying to use the code you gave, modified for my problem.

    md.add_initialized_fem_data("f", mf, density)#Line that I am adding and removing
    md.add_interpolate_transformation_from_expression("V6toS3", m3, m, "f")
    
    I_S3 = gf.asm_generic(mimS3, 0, "Interpolate(f,V6toS3)", -1, md)

The V6 6D volume is the m, mf, md objects, while the objects with S3 (or just 3, such as m3) in them are the ones on the 3D surface, named similarly to how you used them in the code example. f is solved for earlier, and density was computed with

density         = md.variable('f')

, with f being the solved variable from the original 6D FEM method.

Please advise. I think I am doing something incorrectly, but I am not sure what.

Sincerely,
Eric Comstock

6

just use a temporary model for this operation, you do not need to use your main model.

mdtmp = gf.Model("real")
mdtmp.add_initialized_fem_data("somevariable", mfV6, SOMEVARIABLENODALDATA)
mdtmp.add_interpolate_transformation_from_expression("V6toS3", mS3, mV6, "X")

I_S3 = gf.asm_generic(mimS3, 0, "Interpolate(somevariable,V6toS3)", -1, mdtmp)

del mdtmp
7

Hello Kostas,

Thank you!

Should I be using my data extracted with

density         = md.variable('f')

for SOMEVARIABLENODALDATA, or is it supposed to represent something else?

Also, what is "X" supposed to be?

Finally, I have noticed that I_S3 gives a float value (which always seems to be 0.0), not a 3D array. Do you know what could be causing this?

Sincerely,
Eric Comstock

8

you could even write

mdtmp.add_initialized_fem_data("somevariable",
                               md.mesh_variable("f"),
                               md.variable("f"))

this syntax is only for scalars, which means that somevariable in I_S2 should be just one component of the vector, for example “somevariable(1)”, “somevariable(2)”, etc.

X is simply the coordinates vector of the current integration point

9

Hello Kostas,

Ah, I see. What format should I put X in? Should I use numpy format inside the string (e.g. "[1, 1, 2]" for the point x=1, y=1,z=2), or something else?

By the way, I am getting the error message

AttributeError: 'Model' object has no attribute 'mesh_variable'. Did you mean: 'set_variable'?

Would there be performance issues for trying to do this for a 3D grid of points?

Is there any syntax for directly creating a 3D array of the integrals, for each point in space? If so, would that be preferred?

Sincerely,
Eric Comstock

10

no, “X” is just “X”, it is part of the GWFL language

11

Hello Kostas,

I am afraid that I am still confused - I have the code working as

    #### Density integration ####
    
    T_3to6 = np.zeros((6,3))
    T_3to6[3,0] = T_3to6[4,1] = T_3to6[5,2] = 1.
    mS3.transform(T_3to6)
    
    mdtmp = gf.Model("real")
    mdtmp.add_initialized_fem_data("somevariable", mf, density)
    mdtmp.add_interpolate_transformation_from_expression("V6toS3", mS3, m, "X")
    mimS3 = gf.MeshIm(mS3)
    
    I_S3 = gf.asm_generic(mimS3, 0, "Interpolate(somevariable,V6toS3)", -1, mdtmp)
    
    del mdtmp
    logging.debug('Type of I_S3: '+str(type(I_S3)))
    logging.debug('I_S3: '+str(I_S3))

, which runs without errors, but the result of I_S3 is always a float with value 0.0. If "X" is the current integration point over the velocity mesh (mS3), then how do I get the value of the integral as a function of position (either as a getFEM object, or by inputting the position somewhere)?

mS3 is defined previously as a 3D mesh based on the velocity grids.

Sincerely,
Eric Comstock

12

the code looks correct, and should give nonzero if density array has nonzero entries.

To debug this situation, try to replace

I_S3 = gf.asm_generic(mimS3, 0, "Interpolate(somevariable,V6toS3)", -1, mdtmp)

with

I_S3 = gf.asm_generic(mimS3, 0, "Interpolate(Print(somevariable),V6toS3)", -1, mdtmp)

or

I_S3 = gf.asm_generic(mimS3, 0, "Print(Interpolate(somevariable,V6toS3))", -1, mdtmp)

This will print the values of the respective subexpression at each integration point.

13

Hello Kostas,

Thank you! I can confirm that density is nonzero.

I_S3 = gf.asm_generic(mimS3, 0, "Interpolate(Print(somevariable),V6toS3)", -1, mdtmp)

game me the error

RuntimeError: (Getfem::InterfaceError) -- Error in getfem_generic_assembly_tree.cc, line 1679 getfem::GA_TOKEN_TYPE getfem::ga_read_term(getfem::pstring, bgeot::size_type&, getfem::ga_tree&, getfem::ga_macro_dictionary&):
Error in assembly string

and

I_S3 = gf.asm_generic(mimS3, 0, "Print(Interpolate(somevariable,V6toS3))", -1, mdtmp)", -1, mdtmp)

gave me the same I_S3 = 0.0 float that I got before.

Sincerely,
Eric Comstock

14

I see, we forgot to define the integration method, this line is not enough for defining the integration method. Choose one as you did with the other MeshIm objects in your script.

15

Hello Kostas,

Thank you! Apologies for the delay - I spend the day yesterday travelling.

I tried using

mimS3 = gf.MeshIm(m3, gf.Integ('IM_NC(3,' + order + ')'))

for the surface mesh definition - this is the same definition I am using on the 6D mesh (both are simplex-based, which seems to work the best). Note that order is the order of my integration - here it is 1.

Edit: For the 6D mesh I am using 6 dimensions, for this one I am using 3 - just thought I would clarify.

This produced the error


  File "c:\users\ecomstock3\onedrive - georgia institute of technology\documents\getfem simulations\test6d_logfile_profiling.py", line 436, in calc_F
    I_S3 = gf.asm_generic(mimS3, 0, "Interpolate(f,V6toS3)", -1, md)

  File "C:\Users\ecomstock3\AppData\Local\anaconda3\envs\Research_FEM\lib\site-packages\getfem\getfem.py", line 5808, in asm_generic
    return getfem('asm', 'generic', mim, order, expression, region, model, *args)

RuntimeError: (Getfem::InterfaceError) -- Error in getfem_generic_assembly_interpolation.cc, line 663 : 
Incompatibility of meshes in interpolation


logic_error exception caught

Please advise.

Sincerely,
Eric Comstock

16

Hello Kostas,

I it partially works now - I am currently using

    mdtmp = gf.Model("real")
    mdtmp.add_initialized_fem_data("f", mf, density)
    mdtmp.add_interpolate_transformation_from_expression("V6toS3", m3, m, "X")
    
    I_S3 = gf.asm_generic(mimS3, 1, "Interpolate(f,V6toS3)", -1, mdtmp)
    logging.debug('Type of I_S3: ' + str(type(I_S3)))
    logging.debug('Result of density integration vs space: I_S3 = ' + str(I_S3))

for my code (mimS3 and m3 are defined previously), and am not getting any errors. However, the result I am getting is an empty numpy array,

2025-07-03 09:18:17,221 [DEBUG] Type of I_S3: <class 'numpy.ndarray'>
2025-07-03 09:18:17,222 [DEBUG] Result of density integration vs space: I_S3 = []

Do you know what could be causing this?

Sincerely,
Eric Comstock