ISMIP HOM periodic boundary conditions

pf4d

Hello,

I would like to create the periodic boundary conditions appropriate for the ISMIP HOM tests. I see on page 57 of userguide.pdf that I need to create a "pairing" vector for degrees of freedom along certain edges. When querying the data structure, I get

In : md.stressbalance.vertex_pairing

Out: array(, dtype=float64)

I can not find any documentation indicating how to initialize the values; anyone?

mathieumorlighem

Hi Evan,

you can look at examples/ISMIP/CheatyRunme.m. Make sure you use squaremesh (not triangle or bamg). We basically need two columns and each row indicates two vertices that are connected:


   maxX=find(md.mesh.x==max(md.mesh.x));
   minX=find(md.mesh.x==min(md.mesh.x));
   maxY=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=max(md.mesh.x) & md.mesh.x~=min(md.mesh.x));
   minY=find(md.mesh.y==min(md.mesh.y) & md.mesh.x~=max(md.mesh.x) & md.mesh.x~=min(md.mesh.x));
   md.stressbalance.vertex_pairing=;

This is in MATLAB but I am sure you can do something similar with python.

Good luck!

Mathieu

pf4d

Ah ha!

Thank you very much, this is quite clear.

mathieumorlighem

I think you can use something like that:


posx=np.where(md.mesh.x==0.)
posx2=np.where(md.mesh.x==np.max(md.mesh.x))

posy=np.where(np.logical_and.reduce((md.mesh.y==0.,md.mesh.x!=0.,md.mesh.x!=np.max(md.mesh.x)))) #Don't take the same nodes two times
posy2=np.where(np.logical_and.reduce((md.mesh.y==np.max(md.mesh.y),md.mesh.x!=0.,md.mesh.x!=np.max(md.mesh.x))))

md.stressbalance.vertex_pairing=np.vstack((np.vstack((posx+1,posx2+1)).T,np.vstack((posy+1,posy2+1)).T))

pf4d

So the python implementation is fairly straight-forward, I am fairly certain this is correct :

# set periodic lateral boundary conditions :


minX = np.where(md.mesh.x == 0)
maxX = np.where(md.mesh.x == L)

# for y, maxX and minX should be excluded :
minY = np.where(np.logical_and(md.mesh.y == 0,
                               md.mesh.x != L,
                               md.mesh.x != 0))
maxY = np.where(np.logical_and(md.mesh.y == L,
                               md.mesh.x != L,
                               md.mesh.x != 0))

# set the nodes that should be paired together :

md.stressbalance.vertex_pairing = np.array().T

However, there is a strange artifact on the North edge :

pf4d

Ah yes, I see your +1 to the indices, this fixed the issue:

with code:

# set periodic lateral boundary conditions :
minX = np.where(md.mesh.x == 0) + 1 maxX = np.where(md.mesh.x == L) + 1 # for y, maxX and minX should be excluded : minY = np.where(np.logical_and(md.mesh.y == 0, md.mesh.x != L, md.mesh.x != 0)) + 1 maxY = np.where(np.logical_and(md.mesh.y == L, md.mesh.x != L, md.mesh.x != 0)) + 1 # set the nodes that should be paired together :
md.stressbalance.vertex_pairing = np.array().T

I did not see your post because I was formulating a response in this thread for some time.

mathieumorlighem

Yes, md.mesh uses an indexing that starts at +1 (not 0), hence the offset.