Prolong FE function from a coarse to a fine mesh #888

BarbaraV · 2022-04-02T16:16:57Z

BarbaraV
Apr 2, 2022

Hi all,
I am quite new to scikit-fem, so please excuse my (possibly) dumb question.

I am trying to interpret a FE function on a "coarse" mesh like mesh_coarse = skfem.MeshTri().refined(2) also as a FE function on a finer mesh like mesh_fine = skfem.MeshTri().refined(4). How can I accomplish this? I have tried around with Basis.interpolate and Basis.project, but did not succeed.

The background is that I have a FE function on mesh_fine (seen as a reference solution to some problem) and a FE function on mesh_coarse and I would like to compare those two/compute errors between the two.

Thank you very much for your help! I can provide more details on my question if needed.

Answered by kinnala

Apr 3, 2022

In principle, it is possible to use Basis.interpolator between arbitrary meshes.
This is not very efficient though for very fine meshes.

In [1]: from skfem import *                                                     
In [2]: m = MeshTri().refined(2)                                                
In [3]: basis = Basis(m, ElementTriP1())                                        
In [4]: mfine = m.refined(2)                                                    
In [6]: import numpy as np                                                      
In [7]: x = np.sin(5. * m.p[0])                                                 
In [8]: interp = basis.interpolator(x)                                    …

View full answer

sxwang1991 · 2022-04-02T16:31:38Z

sxwang1991
Apr 2, 2022

Did you use Two-grid method?

…

---- 回复的原邮件 ---- | 发件人 | ***@***.***> | | 日期 | 2022年04月03日 00:17 | | 收件人 | ***@***.***> | | 抄送至 | ***@***.***> | | 主题 | [kinnala/scikit-fem] Prolong FE function from a coarse to a fine mesh (Discussion #888) | Hi all, I am quite new to scikit-fem, so please excuse my (possibly) dumb question. I am trying to interpret a FE function on a "coarse" mesh like mesh_coarse = skfem.MeshTri().refined(2) also as a FE function on a finer mesh like mesh_fine = skfem.MeshTri().refined(4). How can I accomplish this? I have tried around with Basis.interpolate and Basis.project, but did not succeed. The background is that I have a FE function on mesh_fine (seen as a reference solution to some problem) and a FE function on mesh_coarse and I would like to compare those two/compute errors between the two. Thank you very much for your help! I can provide more details on my question if needed. — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

0 replies

sxwang1991 · 2022-04-02T16:44:49Z

sxwang1991
Apr 2, 2022

refined_Mesh, X = basis.refinterp(x, 2)

…

---- 回复的原邮件 ---- | 发件人 | ***@***.***> | | 日期 | 2022年04月03日 00:31 | | 收件人 | ***@***.***> | | 抄送至 | | | 主题 | Re: [kinnala/scikit-fem] Prolong FE function from a coarse to a fine mesh (Discussion #888) | Did you use Two-grid method?

---- 回复的原邮件 ---- | 发件人 | ***@***.***> | | 日期 | 2022年04月03日 00:17 | | 收件人 | ***@***.***> | | 抄送至 | ***@***.***> | | 主题 | [kinnala/scikit-fem] Prolong FE function from a coarse to a fine mesh (Discussion #888) | Hi all, I am quite new to scikit-fem, so please excuse my (possibly) dumb question. I am trying to interpret a FE function on a "coarse" mesh like mesh_coarse = skfem.MeshTri().refined(2) also as a FE function on a finer mesh like mesh_fine = skfem.MeshTri().refined(4). How can I accomplish this? I have tried around with Basis.interpolate and Basis.project, but did not succeed. The background is that I have a FE function on mesh_fine (seen as a reference solution to some problem) and a FE function on mesh_coarse and I would like to compare those two/compute errors between the two. Thank you very much for your help! I can provide more details on my question if needed. — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

0 replies

kinnala · 2022-04-03T06:00:43Z

kinnala
Apr 3, 2022
Maintainer

In principle, it is possible to use Basis.interpolator between arbitrary meshes.
This is not very efficient though for very fine meshes.

In [1]: from skfem import *                                                     
In [2]: m = MeshTri().refined(2)                                                
In [3]: basis = Basis(m, ElementTriP1())                                        
In [4]: mfine = m.refined(2)                                                    
In [6]: import numpy as np                                                      
In [7]: x = np.sin(5. * m.p[0])                                                 
In [8]: interp = basis.interpolator(x)                                          
In [10]: xfine = interp(mfine.p)                                                
In [11]: from skfem.visuals.matplotlib import plot3, show                       
In [12]: basisfine = Basis(mfine, ElementTriP1())                               
In [13]: plot3(basisfine, xfine, nrefs=0).show()

The most efficient implementation would create a prolongation matrix during call to MeshTri.refined but this is not currently supported. It shouldn't be hard to implement though. One should just figure out how it fits with the existing functionality.

9 replies

kinnala Apr 7, 2022
Maintainer

There is some room for optimization though:

In [26]: m = MeshTri().refined(8); m.t.shape                                    
Out[26]: (3, 131072)

In [27]: mr = m.refined(); mr.t.shape                                           
Out[27]: (3, 524288)

In [28]: basis = Basis(m, ElementTriP1()); basisr = Basis(mr, ElementTriP1())   

In [29]: %time P = basis.probes(basisr.doflocs)                                 
---------------------------------------------------------------------------
MemoryError                               Traceback (most recent call last)
<timed exec> in <module>

~/src/scikit-fem/skfem/assembly/basis/cell_basis.py in probes(self, x)
    170 
    171         """
--> 172         cells = self.mesh.element_finder(mapping=self.mapping)(*x)
    173         pts = self.mapping.invF(x[:, :, np.newaxis], tind=cells)
    174         phis = np.array(

~/src/scikit-fem/skfem/mesh/mesh_tri_1.py in finder(x, y, _search_all)
    402                 ix = np.arange(nelems, dtype=np.int64)
    403 
--> 404             X = mapping.invF(np.array([x, y])[:, None], ix)
    405             inside = ((X[0] >= 0) *
    406                       (X[1] >= 0) *

~/src/scikit-fem/skfem/mapping/mapping_affine.py in invF(self, x, tind)
    138             invA, b = self.invA[:, :, tind], self.b[:, tind]
    139 
--> 140         y = (x.T - b.T).T
    141 
    142         return np.einsum('ijk,jkl->ikl', invA, y)

MemoryError: Unable to allocate 514. GiB for an array with shape (263169, 131072, 2) and data type float64

kinnala Apr 7, 2022
Maintainer

There should be a faster way of doing this, something along the lines proposed by @sxwang1991 but which takes care of the duplicate nodes introduced by that approach.

xida2020 Jun 4, 2022

Is Basis.interpolator supported for ElementVector(ElementTriP1())?

kinnala Jun 4, 2022
Maintainer

I don't recall testing it. If it doesn't work, you can do basis.split(x) and then use interpolator componentwise.

sxwang1991 Jun 18, 2022

`
from skfem import *
import numpy as np

m = MeshTri().refined(8)
basis = Basis(m, ElementTriP1())
mfine = m.refined(2)

x = np.sin(5. * m.p[0])

mfine_r, xfine_r = basis.refinterp(x, 2)
pts_r = mfine_r.doflocs.T
pts_r_ix = np.unique(pts_r, return_index =True, axis=0)[1]
xfine_r = xfine_r[pts_r_ix,]
pts = mfine.doflocs.T
ix = np.lexsort((pts[:, 1], pts[:, 0]))
ix_inv = np.arange(pts.shape[0])

ix_inv[ix,]=np.arange(pts.shape[0])
xfine_r = xfine_r[ix_inv,]
`

sxwang1991 · 2022-06-18T14:48:40Z

sxwang1991
Jun 18, 2022

I try to use refinterp and it works perfectly and gives the same results as Basis.interpolator.
`
from skfem import *
import numpy as np

m = MeshTri().refined(5)
basis = Basis(m, ElementTriP1())
mfine = m.refined(2)

x = np.sin(5. * m.p[0])
interp = basis.interpolator(x)
xfine = interp(mfine.p)

mfine_r, xfine_r = basis.refinterp(x, 2)

pts_r = mfine_r.doflocs.T
pts_r_ix = np.unique(pts_r, return_index =True, axis=0)[1]
xfine_r = xfine_r[pts_r_ix,]
pts = mfine.doflocs.T
ix = np.lexsort((pts[:, 1], pts[:, 0]))
ix_inv = np.arange(pts.shape[0])

ix_inv[ix,]=np.arange(pts.shape[0])
xfine_r = xfine_r[ix_inv,]

get the same result

print(np.max(xfine-xfine_r))

`

1 reply

sxwang1991 Jun 18, 2022

It returns 0.0.
I haven't yet compared the performance yet.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Prolong FE function from a coarse to a fine mesh #888

{{title}}

Replies: 4 comments 10 replies

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

Select a reply

Prolong FE function from a coarse to a fine mesh #888

BarbaraV Apr 2, 2022

Replies: 4 comments · 10 replies

sxwang1991 Apr 2, 2022

sxwang1991 Apr 2, 2022

kinnala Apr 3, 2022 Maintainer

kinnala Apr 7, 2022 Maintainer

kinnala Apr 7, 2022 Maintainer

xida2020 Jun 4, 2022

kinnala Jun 4, 2022 Maintainer

sxwang1991 Jun 18, 2022

sxwang1991 Jun 18, 2022

get the same result

sxwang1991 Jun 18, 2022

BarbaraV
Apr 2, 2022

Replies: 4 comments 10 replies

sxwang1991
Apr 2, 2022

sxwang1991
Apr 2, 2022

kinnala
Apr 3, 2022
Maintainer

kinnala Apr 7, 2022
Maintainer

kinnala Apr 7, 2022
Maintainer

kinnala Jun 4, 2022
Maintainer

sxwang1991
Jun 18, 2022