# Best practice for a combined integration and interpolation **URL:** \ **Category:** General\ **Created:** [June 26, 2025, 1:27pm UTC](https://getfem.discourse.group/t/best-practice-for-a-combined-integration-and-interpolation/92 "2025-06-26T13:27:14Z")\ **Posts on this page:** 18\ **Page:** 1
**Author:** ![EricAComstock](https://yyz2.discourse-cdn.com/free1/user_avatar/getfem.discourse.group/ericacomstock/32/68_2.png) [@EricAComstock](https://getfem.discourse.group/u/EricAComstock)\ **Post date:** [June 26, 2025, 1:27pm UTC](https://getfem.discourse.group/t/best-practice-for-a-combined-integration-and-interpolation/92/1 "2025-06-26T13:27:14Z")
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 ---
**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\ **Post date:** [June 26, 2025, 5:29pm UTC](https://getfem.discourse.group/t/best-practice-for-a-combined-integration-and-interpolation/92/2 "2025-06-26T17:29:15Z")
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: ```python 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)") ``` ---
**Author:** ![EricAComstock](https://yyz2.discourse-cdn.com/free1/user_avatar/getfem.discourse.group/ericacomstock/32/68_2.png) [@EricAComstock](https://getfem.discourse.group/u/EricAComstock)\ **Post date:** [June 26, 2025, 7:29pm UTC](https://getfem.discourse.group/t/best-practice-for-a-combined-integration-and-interpolation/92/3 "2025-06-26T19:29:32Z")
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 ---
**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\ **Post date:** [June 26, 2025, 9:12pm UTC](https://getfem.discourse.group/t/best-practice-for-a-combined-integration-and-interpolation/92/4 "2025-06-26T21:12:19Z")
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, … ```python 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. ---
**Author:** ![EricAComstock](https://yyz2.discourse-cdn.com/free1/user_avatar/getfem.discourse.group/ericacomstock/32/68_2.png) [@EricAComstock](https://getfem.discourse.group/u/EricAComstock)\ **Post date:** [June 30, 2025, 12:55pm UTC](https://getfem.discourse.group/t/best-practice-for-a-combined-integration-and-interpolation/92/5 "2025-06-30T12:55:22Z")
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 ```auto 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 ```auto 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. ```auto 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 ```auto 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 ---
**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\ **Post date:** [June 30, 2025, 1:54pm UTC](https://getfem.discourse.group/t/best-practice-for-a-combined-integration-and-interpolation/92/6 "2025-06-30T13:54:12Z")
just use a temporary model for this operation, you do not need to use your main model. ```auto 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 ``` ---
**Author:** ![EricAComstock](https://yyz2.discourse-cdn.com/free1/user_avatar/getfem.discourse.group/ericacomstock/32/68_2.png) [@EricAComstock](https://getfem.discourse.group/u/EricAComstock)\ **Post date:** [June 30, 2025, 2:24pm UTC](https://getfem.discourse.group/t/best-practice-for-a-combined-integration-and-interpolation/92/7 "2025-06-30T14:24:48Z")
Hello Kostas, Thank you! Should I be using my data extracted with ```auto 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 ---
**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\ **Post date:** [June 30, 2025, 2:42pm UTC](https://getfem.discourse.group/t/best-practice-for-a-combined-integration-and-interpolation/92/8 "2025-06-30T14:42:45Z")
you could even write ```auto 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 ---
**Author:** ![EricAComstock](https://yyz2.discourse-cdn.com/free1/user_avatar/getfem.discourse.group/ericacomstock/32/68_2.png) [@EricAComstock](https://getfem.discourse.group/u/EricAComstock)\ **Post date:** [June 30, 2025, 3:03pm UTC](https://getfem.discourse.group/t/best-practice-for-a-combined-integration-and-interpolation/92/9 "2025-06-30T15:03:18Z")
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 ```auto 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 ---
**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\ **Post date:** [June 30, 2025, 3:04pm UTC](https://getfem.discourse.group/t/best-practice-for-a-combined-integration-and-interpolation/92/10 "2025-06-30T15:04:40Z")
no, “X” is just “X”, it is part of the GWFL language ---
**Author:** ![EricAComstock](https://yyz2.discourse-cdn.com/free1/user_avatar/getfem.discourse.group/ericacomstock/32/68_2.png) [@EricAComstock](https://getfem.discourse.group/u/EricAComstock)\ **Post date:** [June 30, 2025, 4:25pm UTC](https://getfem.discourse.group/t/best-practice-for-a-combined-integration-and-interpolation/92/11 "2025-06-30T16:25:17Z")
Hello Kostas, I am afraid that I am still confused - I have the code working as ```auto #### 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 ---
**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\ **Post date:** [June 30, 2025, 7:40pm UTC](https://getfem.discourse.group/t/best-practice-for-a-combined-integration-and-interpolation/92/12 "2025-06-30T19:40:39Z")
the code looks correct, and should give nonzero if `density` array has nonzero entries. To debug this situation, try to replace ```auto I_S3 = gf.asm_generic(mimS3, 0, "Interpolate(somevariable,V6toS3)", -1, mdtmp) ``` with ```auto I_S3 = gf.asm_generic(mimS3, 0, "Interpolate(Print(somevariable),V6toS3)", -1, mdtmp) ``` or ```auto 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. ---
**Author:** ![EricAComstock](https://yyz2.discourse-cdn.com/free1/user_avatar/getfem.discourse.group/ericacomstock/32/68_2.png) [@EricAComstock](https://getfem.discourse.group/u/EricAComstock)\ **Post date:** [June 30, 2025, 8:13pm UTC](https://getfem.discourse.group/t/best-practice-for-a-combined-integration-and-interpolation/92/13 "2025-06-30T20:13:06Z")
Hello Kostas, Thank you! I can confirm that `density` is nonzero. ```auto I_S3 = gf.asm_generic(mimS3, 0, "Interpolate(Print(somevariable),V6toS3)", -1, mdtmp) ``` game me the error ```auto 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 ```auto 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 ---
**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\ **Post date:** [July 1, 2025, 5:39am UTC](https://getfem.discourse.group/t/best-practice-for-a-combined-integration-and-interpolation/92/14 "2025-07-01T05:39:34Z")
> [@EricAComstock](#): > > ` mimS3 = gf.MeshIm(mS3)` 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. ---
**Author:** ![EricAComstock](https://yyz2.discourse-cdn.com/free1/user_avatar/getfem.discourse.group/ericacomstock/32/68_2.png) [@EricAComstock](https://getfem.discourse.group/u/EricAComstock)\ **Post date:** [July 2, 2025, 1:07pm UTC](https://getfem.discourse.group/t/best-practice-for-a-combined-integration-and-interpolation/92/15 "2025-07-02T13:07:04Z")
Hello Kostas, Thank you! Apologies for the delay - I spend the day yesterday travelling. I tried using ```auto 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 ```auto 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 ---
**Author:** ![EricAComstock](https://yyz2.discourse-cdn.com/free1/user_avatar/getfem.discourse.group/ericacomstock/32/68_2.png) [@EricAComstock](https://getfem.discourse.group/u/EricAComstock)\ **Post date:** [July 3, 2025, 1:20pm UTC](https://getfem.discourse.group/t/best-practice-for-a-combined-integration-and-interpolation/92/16 "2025-07-03T13:20:04Z")
Hello Kostas, I it partially works now - I am currently using ```auto 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, ```auto 2025-07-03 09:18:17,221 [DEBUG] Type of I_S3: 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 ---
**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\ **Post date:** [July 6, 2025, 5:55pm UTC](https://getfem.discourse.group/t/best-practice-for-a-combined-integration-and-interpolation/92/17 "2025-07-06T17:55:28Z")
the second argument in `gf.asm_generic`: 0 - scalar integral 1 - integrals for all test functions, i.e. residual vector 2 - integrals of combination of 1st and 2nd order test functions, i,e, coupling matrix ---
**Author:** ![EricAComstock](https://yyz2.discourse-cdn.com/free1/user_avatar/getfem.discourse.group/ericacomstock/32/68_2.png) [@EricAComstock](https://getfem.discourse.group/u/EricAComstock)\ **Post date:** [July 7, 2025, 12:30pm UTC](https://getfem.discourse.group/t/best-practice-for-a-combined-integration-and-interpolation/92/18 "2025-07-07T12:30:23Z")
Hello Kostas, Thank you! Sincerely, Eric Comstock