Complement ExtractDofsonBorder for vectorial finite element space

examples/hpddm/Schur-complement-PETSc.edp

@@ -31,7 +31,7 @@ varf vA(u,v) = int2d(Th)(grad(u)' * grad(v)) + on(1, 3, 4, u = 0);
     matrix A = vA(Vh, Vh, tgv = -2);
     Mat dA(A);
     int[int] I;
-    ExtractDofsonBorder(2, Vh, I, 1);
+    ExtractDofsonBorder(2, Vh, I, 1, 1);
     real[int] list(Vh.ndof);
     for [j, ij : I] list[ij] = j + 1;
     real[int, int] S(1, 1);

examples/plugin/Schur-Complement-V4.3.edp

@@ -38,7 +38,7 @@ varf va(u,v) = int2d(Th)(grad(u)'*grad(v))+on(1,3,4,u=0);
  matrix A= va(Vh,Vh,sym=0,tgv=-2);
 if(verbosity>9) cout << A << endl; 
 int[int] I; 
-ExtractDofsonBorder(2,Vh,I,1);
+ExtractDofsonBorder(2,Vh,I,1,1);
 
 real[int,int] B(1,1);
 

examples/plugin/Schur-Complement.edp

@@ -44,7 +44,7 @@ varf va(u,v) = int2d(Th)(grad(u)'*grad(v))+on(1,3,4,u=0);
  matrix A= va(Vh,Vh,sym=0,tgv=-2);
 if(verbosity>9) cout << A << endl; 
 int[int] I; 
-ExtractDofsonBorder(2,Vh,I,1);
+ExtractDofsonBorder(2,Vh,I,1,1);
 
 real[int,int] B(1,1);
 

examples/tutorial/BEM.edp

@@ -41,7 +41,7 @@ Vh<complex> u;
  // computation of the matrice BEM
   // nb of  DoF on border 
  int[int] IdfB2Vh(1); // for numbering IdfB2Vh[i]==i 
- ExtractDofsonBorder(labup,Vh,IdfB2Vh,-1)
+ ExtractDofsonBorder(labup,Vh,IdfB2Vh,-1,1)
  int kdfBEM=IdfB2Vh.n;
  // verif
  if(0)  { Vh X=x;

idp/ExtractDofsonBorder.idp

@@ -10,31 +10,76 @@
    - int[int]  doflabs  (resize array)               (out)  
    - int orient : orientation of the border direct (1) or in direct -1 (in)
 */
-macro ExtractDofsonBorder(labs,Wh,doflabs,orient)
+macro ExtractDofsonBorder(labs,Wh,doflabs,orient,opt)
 {
-mesh Th=Wh.Th;
-fespace VhExtractDofsonBorder(Th,P1);
-VhExtractDofsonBorder abslabs,onlon,atest;
-real lenborder = int1d(Th,labs,qfe=qf1pElump)(1.);
-if( lenborder >0) 
-{
-if(orient)
-solve Pbabciss(abslabs,atest,solver=GMRES) 
-= int2d(Th,qft=qf1pTlump)(abslabs*atest*1e-8/area)
-  +  int1d(Th,labs,qfe=qf1pElump)([dx(abslabs),dy(abslabs)]'*[-N.y,+N.x]*atest)
-  -  int1d(Th,labs,qfe=qf1pElump)(orient*atest);
-  else abslabs= 1.; /* no sorting*/ 
-  real abslabsmax = abs(2*abslabs[].max);
-  varf vlabsneg(u,v) = on(labs,u=abslabs-abslabsmax);
-  if(0)  plot(abslabs,wait=1);
-  real[int] absc=vlabsneg(0,Wh);
-  doflabs.resize(Wh.ndof);
-  doflabs=0:Wh.ndof-1;
-  sort(absc,doflabs);
-  absc= absc? 1:0;
-  doflabs.resize(absc.sum+0.5);
-}
-else 
- doflabs.resize(0);	
+	mesh Th = Wh.Th;
+	real lenborder = int1d(Th,labs,qfe=qf1pElump)(1.);
+	if (lenborder > 0) {
+		if (opt == 1) {
+			fespace VhExtractDofsonBorder(Th,P1);
+			VhExtractDofsonBorder abslabs,onlon,atest;
+			if (orient) {
+				solve Pbabciss(abslabs,atest,solver=GMRES)
+					= int2d(Th,qft=qf1pTlump)(abslabs*atest*1e-8/area)
+					+  int1d(Th,labs,qfe=qf1pElump)([dx(abslabs),dy(abslabs)]'*[-N.y,+N.x]*atest)
+					-  int1d(Th,labs,qfe=qf1pElump)(orient*atest);
+			}
+			else abslabs = 1.; /* no sorting */
+			real abslabsmax = abs(2*abslabs[].max);
+			varf vlabsneg(unused,v) = on(labs,unused=abslabs-abslabsmax);
+			if(0) plot(abslabs,wait=1);
+			real[int] absc=vlabsneg(0,Wh);
+			doflabs.resize(Wh.ndof);
+			doflabs=0:Wh.ndof-1;
+			sort(absc,doflabs);
+			absc= absc? 1:0;
+			doflabs.resize(absc.sum+0.5);
+		} else if (opt == 2) {
+			fespace VhExtractDofsonBorder(Th,[P1,P1]);
+			VhExtractDofsonBorder [abslabs1,abslabs2],[onlon1,onlon2],[atest1,atest2];
+			if (orient) {
+				solve Pbabciss([abslabs1,abslabs2],[atest1,atest2],solver=GMRES)
+					= int2d(Th,qft=qf1pTlump)([abslabs1,abslabs2]'*[atest1,atest2]*1e-8/area)
+					+  int1d(Th,labs,qfe=qf1pElump)([dx(abslabs1),dy(abslabs1)]'*[-N.y,+N.x]*atest1
+							+ [dx(abslabs2),dy(abslabs2)]'*[-N.y,+N.x]*atest2)
+					-  int1d(Th,labs,qfe=qf1pElump)([orient,orient]'*[atest1,atest2]);
+			}
+			else [abslabs1, abslabs2] = [1., 1.]; /* no sorting */
+			real abslabsmax = abs(2*abslabs1[].max);
+			varf vlabsneg([unused1,unused2],[uu,vv])
+				= on(labs,unused1=abslabs1-abslabsmax,unused2=abslabs2-abslabsmax);
+			if(0) plot([abslabs1, abslabs2],wait=1);
+			real[int] absc=vlabsneg(0,Wh);
+			doflabs.resize(Wh.ndof);
+			doflabs=0:Wh.ndof-1;
+			sort(absc,doflabs);
+			absc= absc? 1:0;
+			doflabs.resize(absc.sum+0.5);
+		} else {
+			fespace VhExtractDofsonBorder(Th,[P1,P1,P1]);
+			VhExtractDofsonBorder [abslabs1,abslabs2,abslabs3],[onlon1,onlon2,onlon3],[atest1,atest2,atest3];
+			if (orient) {
+				solve Pbabciss([abslabs1,abslabs2,abslabs3],[atest1,atest2,atest3],solver=GMRES)
+					= int2d(Th,qft=qf1pTlump)([abslabs1,abslabs2,abslabs3]'*[atest1,atest2,atest3]*1e-8/area)
+					+  int1d(Th,labs,qfe=qf1pElump)([dx(abslabs1),dy(abslabs1)]'*[-N.y,+N.x]*atest1
+							+ [dx(abslabs2),dy(abslabs2)]'*[-N.y,+N.x]*atest2
+							+ [dx(abslabs3),dy(abslabs3)]'*[-N.y,+N.x]*atest3)
+					-  int1d(Th,labs,qfe=qf1pElump)([orient,orient,orient]'*[atest1,atest2,atest3]);
+			}
+			else [abslabs1, abslabs2, abslabs3] = [1., 1., 1.]; /* no sorting */
+			real abslabsmax = abs(2*abslabs1[].max);
+			varf vlabsneg([unused1,unused2,unused3],[uu,vv,pp])
+				= on(labs,unused1=abslabs1-abslabsmax,unused2=abslabs2-abslabsmax,unused3=abslabs3-abslabsmax);
+			if(0) plot(abslabs1,wait=1);
+			real[int] absc=vlabsneg(0,Wh);
+			doflabs.resize(Wh.ndof);
+			doflabs=0:Wh.ndof-1;
+			sort(absc,doflabs);
+			absc= absc? 1:0;
+			doflabs.resize(absc.sum+0.5);
+		}
+	} else {
+		doflabs.resize(0);
+	}
 } 
 //