Can physical points tagged using gmsh be called and referred to in eScript?

Asked by Louise Olsen

I have saved Physical Points 1 - 10 in the gmsh mesh using the following gmsh code:

// Borehole radius of 10cm in reservoir surface of 10m x 10m

lx=10.;
ly=10.;

lxh=lx/2.;
lyh=ly/2.;

mesh_size = lx/20.;

P1=newp; Point(P1) = {0,0,0,mesh_size};
P2=newp; Point(P2) = {lx,0,0,mesh_size};
P3=newp; Point(P3) = {lx,ly,0,mesh_size};
P4=newp; Point(P4) = {0,ly,0,mesh_size};

L1=newl; Line(L1) = {P1,P2};
L2=newl; Line(L2) = {P2,P3};
L3=newl; Line(L3) = {P3,P4};
L4=newl; Line(L4) = {P4,P1};

ll1=newll; Line Loop(ll1) = {L1,L2,L3,L4};
// S1=news;
// Plane Surface(S1) = {ll1};

P5=newp; Point(P5)= {lxh,lyh,0,mesh_size};

a_circle=0.1;
ncircle=10;
dtheta=2*Pi/ncircle;

P5=newp; Point(P5)= {lxh,lyh,0,mesh_size};

For i In {1:ncircle:1}
    PCircle[i-1] = newp;
    theta=dtheta*(i-1);
    Point(PCircle[i-1]) = {lxh+a_circle*Cos(theta), lyh+a_circle*Sin(theta),0,mesh_size/100.};

    If (theta<2.*Pi-dtheta/2.)
       Physical Point(i) = {PCircle[i-1]};
       Printf("Physical point %g at angle theta, %g", PCircle[i-1],theta);
    EndIf

    If (i>1)

       C[i-2] = newl;
       Circle(C[i-2]) = {PCircle[i-2],P5,PCircle[i-1]};

    EndIf

EndFor

Cend=newl;
Circle(Cend) = {PCircle[ncircle-1],P5,PCircle[0]};

CLL=newll;

Line Loop(CLL) = {C[],Cend};
// Sborehole=news;
// Plane Surface(Sborehole) = {CLL};

SWhole = news;
Plane Surface(SWhole) = {ll1,CLL};

Physical Surface("SWhole") = {SWhole};

Am I able to visualise/save output for these points only in eScript? Is there a getTaggedValue command similar to setTaggedValue?

This is why I wanted to save Physical points in the gmsh mesh. The following code does not seem to update the physical points:

from esys.escript import *
from esys.finley import ReadGmsh, Rectangle
from esys.escript.linearPDEs import LinearPDE

from esys.weipa import saveVTK

mydomain=ReadGmsh('Borehole2D.msh',2)

x_mesh=Function(mydomain).getX()

r = x_mesh[0]*0.

for i in range(10):
    r.setTaggedValue(i,i+1)

r.expand()
saveVTK("r.vtu", r=r)

Thanks, Louise

Question information

Language:
English Edit question
Status:
Expired
For:
esys-escript Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Launchpad Janitor (janitor) said :
#1

This question was expired because it remained in the 'Open' state without activity for the last 15 days.