Hi! Is there a command in freeFEM where we can list all edges on the mesh? both internal and on the boundary?
Did you found out how to do this?
If you use the “savemesh” command, all the edges are listed up in the mesh file.
but only the boundary edges, not the internal ones
There is “intAllEdges” . In the saved mesh you also get triangles so you could
look at every pair of points for each element in some other code. I wrote
c++ code to load the meshes for my own stuff and it seems to be
pretty simple but its not part of FF and not sure if there is another way
to do what you want.
The idea is simple, use a finite element with one DOF par edge ,
for example P0edge
a exemple
buildedge.edp (377 Bytes)
1 Like
So Eh(int,int) returns an int and the parameters are element and node rahter than
x,y as is the case in Ed(real,real)?
That is good concise syntax but may not be obvious.
Thanks.
You want to have for a space position the edge number,
do can do some thing like:
Eh num;
num[]=0:num[].n-1; // put lhe value of the dof as the dof number. ie : num[][i] == i
so if x,y is on edge the num(x,y) is the dof number but due to roundoff error
use round(num(x,y)) if abs( round(num(x,y)) -num(x,y) ) < 1e-3 .
I’m trying to extract all edges from this exemple. But when I use the comand “extract” it only gives me the boundary edges and not all 5 edges from the mesh Th. Is there a way I can build a meshL with all five edges?
load "msh3"
mesh Th = square(1,1);
int NbTriangles = Th.nt;
plot(Th, wait=1, fill=1);
meshL ThEdge = extract(Th);
ThEdge = trunc(ThEdge, 1, split=1);
plot(ThEdge, wait=1);
fespace Eh(ThEdge, P0);
Eh xx = x, yy=y;
int nDoFl0 = Eh.ndof;
int nDoFl0K = Eh.ndofK;
cout << nDoFl0 << "\n";
Yes the command extract extraterrestrial the boundary mesh.
this is why a given a exemple to build all edges. cf buildedge.edp
and is it possible to build a mesh from this example you provided?
So you want a meshL build for the edge of a mesh, to day not simple way
But why you need that?
Yes, I need a meshL from all the edges of a mesh.
I have a method that is defined on the boundaries so I can build a FE space using edge partition
To day the only is to build a file form my example an read the girl like in example:
3dCurve/extractMeshes.edp
I’ve managed to create a meshL from all the edges, it works perfectly fine when I build a P0 finite element space over this mesh, but when I try with P1, it does not work, the number of DoFs should be 2 per edge, so 10, and it’s giving me 4
Am I doing something wrong?
load "msh3"
int EkOrder = 1;
//------------------Global Mesh---------------------
int globalH = 2^0;
mesh calP = square(globalH,globalH);
int NbTriangles = calP.nt; //number of triangles
int NbVertices = calP.nv; //number of vertices
int NbBndyEdges = calP.nbe; //number of boundary edges
int NbBndy = 3*NbVertices - 3 - NbBndyEdges; //number of edges
int NbIntEdges = NbBndy - NbBndyEdges; //number of internal edges
//-------------------Mesh info---------------------
//-------------------------vector: nodes in elem K
real[int,int] elemNodesGlobal(NbTriangles,3);
for (int i = 0; i < NbTriangles; i++){
for(int j = 0; j < 3; j++)
elemNodesGlobal(i,j) = int(calP[i][j]);
}
//cout << elemNodesGlobal <<"\n";
real[int,int] elemCoordNodesGlobal(NbVertices, 2);
for(int i=0; i<NbVertices; i++){
elemCoordNodesGlobal(i,0) = calP(i).x;
elemCoordNodesGlobal(i,1) = calP(i).y;
}
// cout << elemCoordNodesGlobal << "\n";
// ----------------DoF of V0 space-------------------
fespace V0(calP,P0);
V0 u0;
V0 u0x = x, u0y = y;
int nDoFu0 = V0.ndof; //number of DoF
int nDoFu0K = V0.ndofK; //number of DoF per elem
//function to identify global elements
for(int i = 0; i < NbTriangles; i++){
u0[][i] = i;
}
//----------------edge orientation matrix per elem
int[int,int] edgeOrient(NbTriangles, 3);
for (int k=0; k<NbTriangles;k++){
for(int i=0; i<3;i++){
int i1 = calP[k][(i+1)%3];
int i2 = calP[k][(i+2)%3];
int orientation = i1 < i2 ? 1 : -1;
edgeOrient(k,i) = orientation;
}
}
// cout << "Edge sign" << "\n";
// cout << edgeOrient << "\n";
//------------------------------- local meshes
mesh[int] Th(NbTriangles);
meshL[int] EdgeH(NbTriangles);
//mesh size
real hSubMesh;
int[int] edgeLabel(6);
int[int,int] edgeSignG(NbTriangles, 6);
for(int K = 0; K < NbTriangles; K++){
real coordX0, coordX1, coordX2;
real coordY0, coordY1, coordY2;
for(int i = 0; i< NbVertices; i++){
if(elemNodesGlobal(K,0) == i){
coordX0 = calP(i).x;
coordY0 = calP(i).y;
}
if(elemNodesGlobal(K,1) == i){
coordX1 = calP(i).x;
coordY1 = calP(i).y;
}
if(elemNodesGlobal(K,2) == i){
coordX2 = calP(i).x;
coordY2 = calP(i).y;
}
}
real[int] elemNodesGlobaltoLocal = elemNodesGlobal(K,:);
// cout << elemNodesGlobaltoLocal << "\n";
real[int, int] nodes = [[coordX0, coordY0], [coordX1,coordY1], [coordX2, coordY2]];
border f0(t=0,1){P.x=nodes(2,0)*t + nodes(1,0)*(1-t);P.y=nodes(2,1)*t + nodes(1,1)*(1-t); label=11;};
border f1(t=0,1){P.x=nodes(0,0)*t + nodes(2,0)*(1-t);P.y=nodes(0,1)*t + nodes(2,1)*(1-t); label=22;};
border f2(t=0,1){P.x=nodes(1,0)*t + nodes(0,0)*(1-t);P.y=nodes(1,1)*t + nodes(0,1)*(1-t); label=33;};
mesh ThK = buildmesh(f0(1) + f1(1) + f2(1));
Th[K] = trunc(ThK, abs(u0 -K) < 1e-5, split=2^0);
mesh auxTh = Th[K];
int NbTrianglesLocal = Th[K].nt;
fespace V0k(auxTh, P0);
V0k hk = hTriangle;
hSubMesh = hk[].max;
int sizeVecEdge;
if(EkOrder == 0){
edgeLabel = [11,22,33];
sizeVecEdge = 3;
}
if(EkOrder == 1){
edgeLabel = [11,11,22,22,33,33];
sizeVecEdge = 6;
}
edgeLabel.resize(sizeVecEdge);
for(int i=0; i<edgeLabel.n; i++){
if(edgeLabel(i) == 11){
edgeSignG(K,i) = edgeOrient(K,0);
}
else if(edgeLabel(i) == 22){
edgeSignG(K,i) = edgeOrient(K,1);
}
else if(edgeLabel(i) == 33){
edgeSignG(K,i) = edgeOrient(K,2);
}
}
edgeSignG.resize(NbTriangles,sizeVecEdge);
EdgeH[K] = extract(ThK);
EdgeH[K] = trunc(EdgeH[K], 1, split=2^0);
}
// cout << edgeLabel << "\n";
// cout << edgeSignG << "\n";
meshL Ehg;
for(int K=0; K<NbTriangles; K++){
Ehg = Ehg + EdgeH[K];
}
// plot(Ehg, wait=1);
fespace EdgeSpace(Ehg, P1);
EdgeSpace ex = x, ey = y;
int ndofEh = EdgeSpace.ndof;
int NbEdgesK = EdgeSpace.ndofK;
int NbPartinEdge;
cout << "Number of DoFs: " << ndofEh << "\n";
cout << "Number of DoFs in each edge: " << NbEdgesK << "\n";
You are wrong, for me the for P1, on edges the dof are the vertices not 2 times the edges in this case.
My code more simple I think,
buildedge-BuildMeshL.edp (967 Bytes)