FreeFem 3.5.x
Public Member Functions | Public Attributes | List of all members
fem::FEM Class Reference

this class drives the resolution of the pde using the Finite Element Method. More...

#include <femSolver.hpp>

Public Member Functions

 DECLARE_TYPE (femMesh::femPoint, femPoint)
 Typedefs.
 
 DECLARE_TYPE (femMesh::femTriangle, femTriangle)
 
 DECLARE_TYPE (creal *, cmatptr)
 
 DECLARE_TYPE (float *, matptr)
 
 FEM (femMeshPtr=0, int quadra=0)
 Constructors, destructor and methods. More...
 
 ~FEM ()
 destructor More...
 
float solvePDE (fcts *param, int how)
 solve the PDE More...
 
creal deriv (int m, creal *f, int ksolv, int i)
 
creal convect (creal *f, creal *u1, creal *u2, float dt, int i)
 
creal rhsConvect (creal *f, creal *u1, creal *u2, float dt, int i)
 
creal fctval (creal *f, float x, float y)
 
int getregion (int k)
 
creal gfemuser (creal what, creal *f, int i)
 
creal P1ctoP1 (creal *f, int i)
 
creal prodscalar (creal *f, creal *g)
 
creal ginteg (int, int, int, creal *, creal *, int)
 
creal binteg (int, int, int, creal *, creal *, int)
 
void initvarmat (int how, int flagcomplexe, int N, fcts *param)
 
void assemble (int how, int flagcomplexe, int N, int k, creal *a, creal *b, fcts *param)
 
void solvevarpde (int N, fcts *param, int how)
 

Public Attributes

float * normlx
 
float * normly
 
int N
 

Detailed Description

this class drives the resolution of the pde using the Finite Element Method.

Author
Christophe Prud'homme Chris.nosp@m.toph.nosp@m.e.Pru.nosp@m.dhom.nosp@m.me@an.nosp@m.n.ju.nosp@m.ssieu.nosp@m..fr
See also
femMesh
Version
#$Id: femSolver.hpp,v 1.2 2001/07/12 14:11:57 prudhomm Exp $#

Constructor & Destructor Documentation

◆ FEM()

fem::FEM::FEM ( femMeshPtr  __t = 0,
int  quadra = 0 
)

Constructors, destructor and methods.

default constructor

72 :
73 __mesh( __t ),
74 __quadra( quadra ),
75 bug( 0 ),
76 nhow( 0 ),
77 nhow1( 0 ),
78 nhow2( 0 ),
79 a1c(),
80 rhsQuadra( 0 )
81{
82 int i, k, baux, nquad;
83
84 ns = __mesh->getNumberOfPoints();
85 nt = __mesh->getNumberOfCells();
86 q = __mesh->rp;
87 me = __mesh->tr;
88 ng = __mesh->ng;
89 ngt = __mesh->ngt;
90 bdth = 0;
91 nquad = __quadra ? 3 * nt : ns;
92 for (k = 0; k < nt; k++)
93 for (i = 0; i <= 2; i++)
94 {
95 baux = abss ((me[k][i] - me[k][next[i]]));
96 bdth = (bdth > baux) ? bdth : baux;
97 }
98 a2.destroy();
99 a2.init (nhowmax);
100
101 for( int __i = 0; __i < nhowmax; __i++)
102 {
103 a1c[__i] = 0;
104 }
105
106 area = new float[nt];
107 normlx = new float[nquad];
108 normly = new float[nquad];
109 for (i = 0;i < nquad;i++)
110 {
111 normlx[i] = 0.F;
112 normly[i] = 0.F;
113 }
114 nhow1 = 0;
115 nhow2 = 0;
116 connectiv ();
117 flag.fem = 1;
118 doedge ();
119 buildarea();
120}

◆ ~FEM()

fem::FEM::~FEM ( )

destructor

123{
124 int i;
125
126 if (flag.fem)
127 {
128 for (i = 0; i < nhow2; i++)
129 a2[i].destroy ();
130 for (i = 0; i < nhow1; i++)
131 {
132 if ( a1[i] == 0 )
133 {
134 delete [] a1[i];
135 a1[i]=0;
136 }
137
138 }
139 for (i = 0; i < nhow1; i++)
140 {
141 if ( a1c[i] == 0 )
142 {
143 delete [] a1c[i];
144 a1c[i]=0;
145 }
146 }
147
148 a2.destroy ();
149 // a1.destroy();
150 delete [] area; area = NULL;
151 delete [] normlx; normlx = NULL;
152 delete [] normly; normly = NULL;
153 delete [] triauptr; triauptr = NULL;
154 delete [] triaunb; triaunb = NULL;
155 delete [] lowV; lowV = NULL;
156 delete [] highV; highV = NULL;
157 delete [] Tleft; Tleft = NULL;
158 delete [] Tright; Tright = NULL;
159 delete [] edgeT; edgeT = NULL;
160 delete [] listHead; listHead = NULL;
161
162 flag.fem = 0;
163 nhow = 0;
164 nhow1 = 0;
165 nhow2 = 0;
166 }
167}

Member Function Documentation

◆ assemble()

void fem::FEM::assemble ( int  how,
int  flagcomplexe,
int  N,
int  k,
creal a,
creal b,
fcts param 
)
1362{
1363 int i,j,in,jn, jglob;
1364 long nsl = ns;
1365 switch (N)
1366 {
1367 case 1:
1368 if (!flagcomplexe)
1369 for(j=0;j<3;j++)
1370 { jglob = me[k][j];
1371 param->sol1[jglob] -= realpart(b[j]);
1372 if(how>0)
1373 for(i=0;i<3;i++)
1374 a1[how-1][nsl*(me[k][i]-jglob+bdth)+jglob]
1375 += realpart(a[i*3+j]);
1376 }
1377 break;
1378 case 2:
1379 for(j=0;j<3;j++)
1380 { jglob = me[k][j];
1381 for(jn=0;jn<N;jn++)
1382 { param->sol2[jglob][jn] -= realpart(b[3*jn+j]);
1383 if(how>0)
1384 for(i=0;i<3;i++)
1385 for(in=0;in<N;in++)
1386 a2[how-1][nsl*(me[k][i]-jglob+bdth)+jglob][in*2+jn]
1387 += realpart(a[18*in+9*jn+i*3+j]);
1388 }
1389 }
1390 break;
1391 }
1392
1393}

◆ binteg()

creal fem::FEM::binteg ( int  ref1,
int  ref2,
int  ref3,
creal f,
creal g,
int  ksolv 
)
1135{ // computes a boundary integral, global if ksolv=1, on one T if ksolv>1
1136 long k;
1137 creal integrale = 0.;
1138
1139 if(ksolv>1)
1140 integrale = binteg_t(ksolv-2, ref1, ref2, ref3, f, g);
1141 else
1142 for (k = 0;k < nt;k++)
1143 integrale += binteg_t(k, ref1, ref2, ref3, f, g);
1144
1145 return integrale;
1146}

◆ convect()

creal fem::FEM::convect ( creal f,
creal u1,
creal u2,
float  dt,
int  i 
)
383 {
384 int j, j1, k, k1, iloc, err = 0;
385 float dt1, xl0, xl1, xl2, x, y, qcx, qcy;
386
387 if (!__quadra)
388 {
389 k = listHead[i1];
390 j = i1;
391 if (k <= 0)
392 return f[i1];
393 else
394 {
395 qcx = (q[me[k][0]][0] + q[me[k][1]][0] + q[me[k][2]][0]) / 3;
396 x = qcx + 0.99F * (q[j][0] - qcx);
397 qcy = (q[me[k][0]][1] + q[me[k][1]][1] + q[me[k][2]][1]) / 3;
398 y = qcy + 0.99F * (q[j][1] - qcy);
399 dt1 = dt;
400 k1 = k;
401 if (0 < xtoX (u1, u2, &dt1, &x, &y, &k1))
402 err++;
403 if (barycoor (x, y, k1, &xl0, &xl1, &xl2))
404 err += 3;
405 return f[me[k1][0]] * xl0 + f[me[k1][1]] * xl1 + f[me[k1][2]] * xl2;
406 }
407 }
408 else //__quadra==1
409
410 {
411 if (i1 == 3 * nt - 1)
412 {
413 creal aux = gconvect[i1];
414 delete [] gconvect; gconvect = NULL;
415 return aux;
416 }
417 else if (i1 > 0)
418 return gconvect[i1];
419 else // i1==0
420
421 {
422 creal gloc[3];
423
424 gconvect = new creal[3*nt];
425 for (k = 0; k < nt; k++)
426 {
427 qcx = (q[me[k][0]][0] + q[me[k][1]][0] + q[me[k][2]][0]) / 3;
428 qcy = (q[me[k][0]][1] + q[me[k][1]][1] + q[me[k][2]][1]) / 3;
429 for (iloc = 0; iloc < 3; iloc++)
430 {
431 j = me[k][iloc];
432 j1 = me[k][next[iloc]];
433 x = qcx + 0.999F * ((q[j][0] + q[j1][0]) / 2 - qcx);
434 y = qcy + 0.999F * ((q[j][1] + q[j1][1]) / 2 - qcy);
435 dt1 = dt;
436 k1 = k;
437 if (0 < xtoX (u1, u2, &dt1, &x, &y, &k1))
438 err++;
439 if (barycoor (x, y, k1, &xl0, &xl1, &xl2))
440 err += 3;
441 gloc[next[iloc + 1]] = f[3 * k1] * xl0 + f[3 * k1 + 1] * xl1 + f[3 * k1 + 2] * xl2;
442 }
443 for (iloc = 0; iloc < 3; iloc++)
444 gconvect[3 * k + iloc] = gloc[next[iloc]] + gloc[next[iloc + 1]] - gloc[iloc];
445 }
446 return gconvect[0];
447 }
448 }
449}

◆ deriv()

creal fem::FEM::deriv ( int  m,
creal f,
int  ksolv,
int  i 
)
1232{
1233 long j, k, l;
1234 creal dfxk = 0.F;
1235 creal sigma = 0.F;
1236
1237 if(ksolv>1)
1238 {
1239 k = ksolv - 2;
1240 for (j = 0; j <= 2; j++)
1241 if (m == 0)
1242 dfxk += dwx (j, k) * f[me[k][j]];
1243 else
1244 dfxk += dwy (j, k) * f[me[k][j]];
1245 return dfxk;
1246 }
1247 if (__quadra)
1248 {
1249 k = i / 3;
1250 for (j = 0; j <= 2; j++)
1251 if (m == 0)
1252 dfxk += dwx (j, k) * f[3 * k + j];
1253 else
1254 dfxk += dwy (j, k) * f[3 * k + j];
1255 return dfxk;
1256 }
1257 else
1258 {
1259 for (l = triauptr[i]; l <= triauptr[i + 1] - 1; l++)
1260 sigma += area[triaunb[l]];
1261 for (l = triauptr[i]; l <= triauptr[i + 1] - 1; l++)
1262 {
1263 k = triaunb[l];
1264 for (j = 0; j <= 2; j++)
1265 if (m == 0)
1266 dfxk += dwxa (j, k) * f[me[k][j]];
1267 else
1268 dfxk += dwya (j, k) * f[me[k][j]];
1269 }
1270 return dfxk / sigma;
1271 }
1272}

◆ fctval()

creal fem::FEM::fctval ( creal f,
float  x,
float  y 
)
573{
574 int j, k, kmin = -1, count = 0, notok = 0;
575 float xc, yc, dist, xl0, xl1, xl2, distmin = (float)1e10;
576 float x3 = 3 * x;
577 float y3 = 3 * y;
578 double mu, nu, xcx, ycy;
579
580 for (k = 0; k < nt; k++)
581 {
582 xc = (q[me[k][0]][0] + q[me[k][1]][0] + q[me[k][2]][0]);
583 yc = (q[me[k][0]][1] + q[me[k][1]][1] + q[me[k][2]][1]);
584 dist = (float)(fabs (xc - x3) + fabs (yc - y3));
585 if (dist < distmin)
586 {
587 kmin = k;
588 distmin = dist;
589 }
590 }
591 xc = (q[me[kmin][0]][0] + q[me[kmin][1]][0] + q[me[kmin][2]][0]) / 3;
592 yc = (q[me[kmin][0]][1] + q[me[kmin][1]][1] + q[me[kmin][2]][1]) / 3;
593 while ((count++ < 20) && (notok = barycoor (x, y, kmin, &xl0, &xl1, &xl2), notok))
594 {
595 xcx = xc - x;
596 ycy = yc - y;
597 j = Tconvect (kmin, xcx, ycy, xc, yc, &mu, &nu);
598 if (j >= 0)
599 {
600 kmin = (kmin == Tleft[j = edgeT[kmin][next[next[j]]]]) ? Tright[j] : Tleft[j];
601 xc += (float)(mu = mmax (mu, -1.0)) * (xc - x);
602 yc += (float)mu * (yc - y);
603 }
604 else
605 count = 200;
606 }
607 if (!notok)
608 {
609 if (__quadra)
610 return f[3 * kmin] * xl0 + f[3 * kmin + 1] * xl1 + f[3 * kmin + 2] * xl2;
611 else
612 return f[me[kmin][0]] * xl0 + f[me[kmin][1]] * xl1 + f[me[kmin][2]] * xl2;
613 }
614 else
615 return (float)2e30;
616}

◆ getregion()

int fem::FEM::getregion ( int  k)
Returns
the ngt of a femTriangle to which belongs vertex k
1277{
1278 return __mesh->ngt[listHead[i]];
1279
1280}

◆ gfemuser()

creal fem::FEM::gfemuser ( creal  what,
creal f,
int  i 
)
1284{ // solves A_0 f = A_1 f
1285
1286 int i, k;
1287 float s;
1288
1289 if (j == 0)
1290 {
1291 float *ff;
1292 ff = new float [ns];
1293
1294 for (i = 0; i < ns; i++)
1295 ff[i] = realpart (f[i]);
1296 for (i = 0; i < ns; i++)
1297 {
1298 s = ff[i];
1299 for (k = i + 1; k <= mmin (ns - 1, i + bdth); k++)
1300 s += a1[1][ns * (i - k + bdth) + k] * ff[k];
1301 ff[i] = s;
1302 }
1303 for (i = ns - 1; i >= 0; i--)
1304 {
1305 s = 0.F;
1306 for (k = mmax (i - bdth, 0); k <= i; k++)
1307 s += a1[1][ns * (i - k + bdth) + k] * ff[k];
1308 ff[i] = s;
1309 }
1310 gaussband (a1[0], ff, ns, bdth, 0, 1.F / penal);
1311 for (i = 0; i < ns; i++)
1312 f[i] = ff[i];
1313 delete [] ff; ff = NULL;
1314 }
1315 return f[j];
1316}

◆ ginteg()

creal fem::FEM::ginteg ( int  ref1,
int  ref2,
int  ref3,
creal f,
creal g,
int  ksolv 
)
1182{// computes a volume integral, global if ksolv=1, on one T if ksolv>1
1183 long k;
1184 creal integrale = 0;
1185
1186 if(ksolv>1)
1187 {
1188 k=ksolv-2;
1189 if (ref1 == 0)
1190 integrale += ginteg_t(k, f, g);
1191 else if (ref2 == 0)
1192 {
1193 if (ngt[k] == ref1)
1194 integrale += ginteg_t(k, f, g);
1195 }
1196 else if (ref3 == 0)
1197 {
1198 if (ngt[k] == ref1 || ngt[k] == ref2)
1199 integrale += ginteg_t(k, f, g);
1200 }
1201 else
1202 if (ngt[k] == ref1 || ngt[k] == ref2 || ngt[k] == ref3)
1203 integrale += ginteg_t(k, f, g);
1204 return integrale;
1205 }
1206 if (ref1 == 0)
1207 for (k = 0;k < nt;k++)
1208 integrale += ginteg_t(k, f, g);
1209 else if (ref2 == 0)
1210 for (k = 0;k < nt;k++)
1211 {
1212 if (ngt[k] == ref1)
1213 integrale += ginteg_t(k, f, g);
1214 }
1215 else if (ref3 == 0)
1216 {
1217 for (k = 0;k < nt;k++)
1218 if (ngt[k] == ref1 || ngt[k] == ref2)
1219 integrale += ginteg_t(k, f, g);
1220 }
1221 else
1222 for (k = 0;k < nt;k++)
1223 if (ngt[k] == ref1 || ngt[k] == ref2 || ngt[k] == ref3)
1224 integrale += ginteg_t(k, f, g);
1225
1226 return integrale;
1227}

◆ initvarmat()

void fem::FEM::initvarmat ( int  how,
int  flagcomplexe,
int  N,
fcts param 
)
1320{
1321 long nsl = ((long) ns) * (2 * ((long) bdth) + 1);
1322 long i;
1323 int factorize = (how>0);
1324
1325 if (how > nhowmax)
1326 erreur ("Too many linear systems");
1327 if((N==1)&&(how >nhow1+1)) erreur ("Linear systems number must be created in the natural order");
1328 if((N==2)&&(how >nhow2+1)) erreur ("Linear systems number must be created in the natural order");
1329 if (how < 0)
1330 {
1331 how = -how;
1332 if (((how > nhow1) && (N == 1)) || ((how > nhow2) && (N == 2)))
1333 {
1334 sprintf (errbuf, "solve(..,'-%d') refers to an inexistant system", how);
1335 erreur (errbuf);
1336 }
1337 }
1338 switch (N)
1339 {
1340 case 1:
1341 if (flag.complexe)
1342 { if (how > nhow1) { a1c[nhow1] = new creal[nsl]; nhow1++; }
1343 if(factorize) for (i = 0; i < nsl; i++) a1c[how-1][i] = 0.F;
1344 }
1345 else
1346 { if (how > nhow1) { a1[nhow1] = new float[nsl]; nhow1++; }
1347 if(factorize) for (i = 0; i < nsl; i++) a1[how-1][i] = 0.F;
1348 for (i = 0; i < ns; i++) param->sol1[i] = 0.F;
1349 }
1350 break;
1351 case 2:
1352 if (how > nhow2) { a2[nhow2].init(nsl); nhow2++;}
1353 if(factorize) for (i = 0; i < nsl; i++)
1354 a2[how-1][i] = 0.F;
1355 for (i = 0; i < ns; i++) param->sol2[i] = 0.F;
1356 break;
1357 }
1358}

◆ P1ctoP1()

creal fem::FEM::P1ctoP1 ( creal f,
int  i 
)
1064{
1065 long j, k, l;
1066 creal meanf = 0.F;
1067
1068 for (l = triauptr[i]; l <= triauptr[i + 1] - 1; l++)
1069 {
1070 k = triaunb[l];
1071 for (j = 0; j <= 2; j++)
1072 if (i == me[k][j])
1073 meanf += f[3 * k + j];
1074 }
1075 return meanf / (float)(triauptr[i + 1] - triauptr[i]);
1076}

◆ prodscalar()

creal fem::FEM::prodscalar ( creal f,
creal g 
)
548{
549 creal prod = 0.F;
550 int i1, i2;
551
552 for (int k = 0; k < nt; k++)
553 for (int iloc = 0; iloc < 3; iloc++)
554 {
555 if (__quadra)
556 {
557 i1 = 3 * k + iloc;
558 i2 = 3 * k + next[iloc];
559 }
560 else
561 {
562 i1 = me[k][iloc];
563 i2 = me[k][next[iloc]];
564 }
565 prod += (f[i1] + f[i2]) * (g[i1].conjug() + g[i2].conjug()) * area[k];
566 }
567 return prod / 12.F;
568}

◆ rhsConvect()

creal fem::FEM::rhsConvect ( creal f,
creal u1,
creal u2,
float  dt,
int  i 
)
454 {
455 int j, j1, k, k1, ke, iloc, err = 0;
456 float dt1, xl0, xl1, xl2, x, y, qcx, qcy, airne;
457 creal gconvectke;
458
459 if (!__quadra)
460 {
461 if (i1 == ns)
462 {
463 creal aux = gconvect[i1];
464
465 delete [] gconvect;gconvect = NULL;
466 return aux;
467 }
468 else if (i1 > 0)
469 return gconvect[i1];
470 else // i1==0
471 {
472 gconvect = new creal[ns];
473 rhsQuadra = 1;
474 for (ke = 0; ke < ne; ke++)
475 {
476 airne = 0.F;
477 if (k1 = Tleft[ke], k1 >= 0)
478 airne += area[k1];
479 if (k = Tright[ke], k >= 0)
480 airne += area[k];
481 else
482 k = k1;
483 qcx = (q[me[k][0]][0] + q[me[k][1]][0] + q[me[k][2]][0]) / 3;
484 qcy = (q[me[k][0]][1] + q[me[k][1]][1] + q[me[k][2]][1]) / 3;
485 j = lowV[ke];
486 j1 = highV[ke];
487 x = qcx + 0.999F * ((q[j][0] + q[j1][0]) / 2 - qcx);
488 y = qcy + 0.999F * ((q[j][1] + q[j1][1]) / 2 - qcy);
489 dt1 = dt;
490 k1 = k;
491 if (0 < xtoX (u1, u2, &dt1, &x, &y, &k1))
492 err++;
493 if (barycoor (x, y, k1, &xl0, &xl1, &xl2))
494 err += 3;
495 gconvectke = f[me[k1][0]] * xl0 + f[me[k1][1]] * xl1 + f[me[k1][2]] * xl2;
496 gconvect[lowV[ke]] += airne * gconvectke / 6.F;
497 gconvect[highV[ke]] += airne * gconvectke / 6.F;
498 }
499 return gconvect[0];
500 }
501 }
502 else //quadra==1
503
504 {
505 if (i1 == 3 * nt - 1)
506 {
507 creal aux = gconvect[i1];
508
509 delete [] gconvect;gconvect = NULL;
510 return aux;
511 }
512 else if (i1 > 0)
513 return gconvect[i1];
514 else // i1==0
515
516 {
517 creal gloc[3];
518
519 gconvect = new creal[ns];
520 for (k = 0; k < nt; k++)
521 {
522 qcx = (q[me[k][0]][0] + q[me[k][1]][0] + q[me[k][2]][0]) / 3;
523 qcy = (q[me[k][0]][1] + q[me[k][1]][1] + q[me[k][2]][1]) / 3;
524 for (iloc = 0; iloc < 3; iloc++)
525 {
526 j = me[k][iloc];
527 j1 = me[k][next[iloc]];
528 x = qcx + 0.999F * ((q[j][0] + q[j1][0]) / 2 - qcx);
529 y = qcy + 0.999F * ((q[j][1] + q[j1][1]) / 2 - qcy);
530 dt1 = dt;
531 k1 = k;
532 if (0 < xtoX (u1, u2, &dt1, &x, &y, &k1))
533 err++;
534 if (barycoor (x, y, k1, &xl0, &xl1, &xl2))
535 err += 3;
536 gloc[next[iloc + 1]] = f[3 * k1] * xl0 + f[3 * k1 + 1] * xl1 + f[3 * k1 + 2] * xl2;
537 }
538 for (iloc = 0; iloc < 3; iloc++)
539 gconvect[3 * k + iloc] = gloc[next[iloc]] + gloc[next[iloc + 1]] - gloc[iloc];
540 }
541 return gconvect[0];
542 }
543 }
544}

◆ solvePDE()

float fem::FEM::solvePDE ( fcts param,
int  how 
)

solve the PDE

Parameters
paramcontain all the possible data for computation
howdefines if the P1 quadrature
1001{
1002 long nsl = ((long) ns) * (2 * ((long) bdth) + 1);
1003 int factorize = 1;
1004
1005 if (how > nhowmax)
1006 erreur ("Too many linear systems");
1007 if (how < 0)
1008 {
1009 factorize = 0;
1010 how = -how;
1011 if (((how > nhow1) && (N == 1)) || ((how > nhow2) && (N == 2)))
1012 {
1013 sprintf (errbuf, "solve(..,'-%d') refers to an inexistant system", how);
1014 erreur (errbuf);
1015 }
1016 }
1017 if (((how > nhow1) && (N == 1)) || ((how > nhow2) && (N == 2)))
1018 {
1019 switch (N)
1020 {
1021 case 1:
1022 if (flag.complexe)
1023 a1c[nhow1++] = new creal[nsl];
1024 else
1025 a1[nhow1++] = new float[nsl];
1026
1027 break;
1028 case 2:
1029 a2[nhow2++].init (nsl);
1030 break;
1031 }
1032 }
1033 if (flag.complexe)
1034 {
1035 if (N == 1)
1036 return pdeian (a1c[how - 1], param->sol1c, param->f1c, param->g1c, param->p1c, param->b1c,
1037 param->nuxx1c, param->nuxy1c, param->nuyx1c, param->nuyy1c, param->a11c,
1038 param->a21c, param->c1c, factorize);
1039 else if (N == 2)
1040 return pdeian (a2[how - 1], param->sol2, param->f2, param->g2, param->p2, param->b2,
1041 param->nuxx2, param->nuxy2, param->nuyx2, param->nuyy2, param->a12,
1042 param->a22, param->c2, factorize);
1043 else
1044 return -1.F;
1045
1046 }
1047 else
1048 {
1049 if (N == 1)
1050 return pdeian (a1[how - 1], param->sol1, param->f1, param->g1, param->p1, param->b1,
1051 param->nuxx1, param->nuxy1, param->nuyx1, param->nuyy1, param->a11,
1052 param->a21, param->c1, factorize);
1053 else if (N == 2)
1054 return pdeian (a2[how - 1], param->sol2, param->f2, param->g2, param->p2, param->b2,
1055 param->nuxx2, param->nuxy2, param->nuyx2, param->nuyy2, param->a12,
1056 param->a22, param->c2, factorize);
1057 else
1058 return -1.F;
1059 }
1060}

◆ solvevarpde()

void fem::FEM::solvevarpde ( int  N,
fcts param,
int  how 
)
1397{
1398 long nsl = ns;
1399 int i, factorize = (how>0);
1400 long mekj, nquad = ns;//__quadra ? 3 * nt : ns;
1401 if(how<0)how = -how;
1402 if(N==1)
1403 {
1404 for (i = 0; i < nquad; i++)
1405 if (norm2 (param->p1[i]) != 0)
1406 {
1407 /* if (__quadra)
1408 {
1409 k = i / 3;
1410 j = i - 3 * k;
1411 mekj = me[k][j];
1412 }
1413 else*/
1414 mekj = i;
1415 param->sol1[mekj] += param->p1[i] * penal;
1416 if (factorize)
1417 {
1418 a1[how-1][nsl * bdth + mekj] += id(param->p1[i]) * penal;
1419 }
1420 gaussband (a1[how-1], param->sol1, nsl, bdth, factorize, (float)(1.0 / penal));
1421 }
1422 else if(N==2)
1423 {
1424 for (i = 0; i < nquad; i++)
1425 if (norm2 (param->p2[i]) != 0)
1426 {
1427 param->sol2[i] += param->p2[i] * penal;
1428 if (factorize)
1429 a2[how-1][nsl * bdth + i] += fem::id(param->p2[i]) * penal;
1430 }
1431 gaussband (a2[how-1], param->sol2, nsl, bdth, factorize, (float)(1.0 / penal));
1432 }
1433 }
1434}

The documentation for this class was generated from the following files:

This is the FreeFEM reference manual
Provided by The KFEM project