Error in 1D integration

General Discussion
1

coarse_grid.mesh (266 Bytes)
Plane_1.mesh (8.7 KB)

Hello all,

I am running this code:

verbosity=0;
load "Element_PkEdge";
load "msh3";

meshS TH=readmeshS("coarse_grid.mesh");
fespace Tri(TH,P0); // P0 on coarse mesh
Tri ChiK;

ChiK[][0]=1; // P0 function used to mark the element i
meshS THK=trunc(TH, ChiK>0.1, split=1); // a mesh made of only one triangle;

int[int] lab = labels(THK);
int[int] NewLab = [lab[0],10,lab[1],20,lab[2],30];
THK = change(THK, label=NewLab);
plot(THK, wait=1);
meshS ThK=readmeshS("Plane_1.mesh");

/*********************************/
meshL THL= extract(THK); // border mesh
fespace Ehc(THL,P0); //  Constant fonction of edge
Ehc ul = region+0.5; // get the label on the each boundary edge
ThK = change(ThK, flabel= dist(THL) < 1e-6 ? ul : 40 ) ;
/**************************************/

fespace VHKP0edge(THK,P0edge); //useful to locate the edges of the coarse element
fespace VhK(ThK, [P2, P2]); //FE space to solve the local problem

varf vA([u,v],[uu, vv])=int2d(ThK)((dx(u)*dx(uu)+dy(u)*dy(uu)+dx(v)*dx(vv)+dy(v)*dy(vv))) + on(40,u=0, v=0);
matrix A=vA(VhK,VhK);

func w1 = 1*[1,0];
varf vL([u], [uu,vv]) = int1d(ThK,10,20,30)(u*w1'*[uu,vv]);
matrix L1=vL(VHKP0edge, VhK);   //Line 33 The problem is here

But I get the following error:

  current line = 33
Assertion fail : (0)
	line :2764, in file problem.cpp
Assertion fail : (0)
	line :2764, in file problem.cpp
 err code 6 ,  mpirank 0

The problem comes from matrix L1=vL(VHKP0edge, VhK).

It’s seems that I can not do 1D integration on this mesh. I don’t know what is the problem.

Could someone help me to understand what is the problem ? I have attached the two meshes used, to do some tests.

Thank you in advance,

Best regards,

Loïc,

4

According to the documentation: int1D() does not support meshS Type.

I have solved my problem by converting my surface meshS into a 2D mesh.

Loïc,

5

Yes, I miss a piece of code, the simple way is the do in 2d, I will add in next version this peace of code .

THE CORRECTION /

MBP-M1-FH:ff-10 hecht$ git diff
diff --git a/src/fflib/problem.cpp b/src/fflib/problem.cpp
index 3d14aacd..b701f528 100755
--- a/src/fflib/problem.cpp
+++ b/src/fflib/problem.cpp
@@ -2760,13 +2760,14 @@ template<class R>
             }
         }
         else // int on edge ie
-        //{
-            ffassert(0);
-            /*R2 PA,PB;
+        {
+          //  ffassert(0);
+            R2 PA,PB;
             PA=TriangleHat[VerticesOfTriangularEdge[ie][0]];
             PB=TriangleHat[VerticesOfTriangularEdge[ie][1]];
             R3 E=T.Edge(ie);
             double le = sqrt((E,E));
+        
             for (npi=0;npi<FIb.n;npi++) // loop on the integration point
             {
                 QuadratureFormular1dPoint pi( FIb[npi]);
@@ -2865,7 +2866,7 @@ template<class R>
                             }
                     }
                 }}
-            }*/
+            }
 
 
         *MeshPointStack(stack) = mp;
MBP-M1-FH:ff-10 hecht$