1 /****************************************************************************/
2 /* This file is part of FreeFEM.                                            */
3 /*                                                                          */
4 /* FreeFEM is free software: you can redistribute it and/or modify          */
5 /* it under the terms of the GNU Lesser General Public License as           */
6 /* published by the Free Software Foundation, either version 3 of           */
7 /* the License, or (at your option) any later version.                      */
8 /*                                                                          */
9 /* FreeFEM is distributed in the hope that it will be useful,               */
10 /* but WITHOUT ANY WARRANTY; without even the implied warranty of           */
11 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the            */
12 /* GNU Lesser General Public License for more details.                      */
13 /*                                                                          */
14 /* You should have received a copy of the GNU Lesser General Public License */
15 /* along with FreeFEM. If not, see <http://www.gnu.org/licenses/>.          */
16 /****************************************************************************/
17 // SUMMARY : ...
18 // LICENSE : LGPLv3
19 // ORG     : LJLL Universite Pierre et Marie Curie, Paris, FRANCE
20 // AUTHORS : ...
21 // E-MAIL  : ...
22 
23 // Example C++ function "myfunction", dynamically loaded into "load.edp"
24 
25 /* clang-format off */
26 //ff-c++-LIBRARY-dep:
27 //ff-c++-cpp-dep:
28 /* clang-format on */
29 
30 #include "ff++.hpp"
31 
32 class MatrixUpWind0 : public E_F0 {
33  public:
34   typedef Matrice_Creuse< R > *Result;
35   Expression emat, expTh, expc, expu1, expu2;
MatrixUpWind0(const basicAC_F0 & args)36   MatrixUpWind0(const basicAC_F0 &args) {
37     args.SetNameParam( );
38     emat = args[0];
39     expTh = to< pmesh >(args[1]);
40     expc = CastTo< double >(args[2]);
41     const E_Array *a = dynamic_cast< const E_Array * >((Expression)args[3]);
42     if (!a || a->size( ) != 2) {
43       CompileError("syntax:  MatrixUpWind0(Th,rhi,[u1,u2])");
44     }
45 
46     int err = 0;
47     expu1 = CastTo< double >((*a)[0]);
48     expu2 = CastTo< double >((*a)[1]);
49   }
50 
~MatrixUpWind0()51   ~MatrixUpWind0( ) {}
52 
typeargs()53   static ArrayOfaType typeargs( ) {
54     return ArrayOfaType(atype< Matrice_Creuse< R > * >( ), atype< pmesh >( ), atype< double >( ),
55                         atype< E_Array >( ));
56   }
57 
f(const basicAC_F0 & args)58   static E_F0 *f(const basicAC_F0 &args) { return new MatrixUpWind0(args); }
59 
60   AnyType operator( )(Stack s) const;
61 };
62 
63 class MatrixUpWind3 : public E_F0 {
64  public:
65   typedef Matrice_Creuse< R > *Result;
66   Expression emat, expTh, expc, expu1, expu2, expu3;
MatrixUpWind3(const basicAC_F0 & args)67   MatrixUpWind3(const basicAC_F0 &args) {
68     args.SetNameParam( );
69     emat = args[0];
70     expTh = to< pmesh3 >(args[1]);
71     expc = CastTo< double >(args[2]);
72     const E_Array *a = dynamic_cast< const E_Array * >((Expression)args[3]);
73     if (a == NULL) printf("Dynamic cast failed\n");
74     if (a->size( ) != 3) {
75       CompileError("syntax:  MatrixUpWind0(Th,rhi,[u1,u2])");
76     }
77 
78     int err = 0;
79     expu1 = CastTo< double >((*a)[0]);
80     expu2 = CastTo< double >((*a)[1]);
81     expu3 = CastTo< double >((*a)[2]);
82   }
83 
~MatrixUpWind3()84   ~MatrixUpWind3( ) {}
85 
typeargs()86   static ArrayOfaType typeargs( ) {
87     return ArrayOfaType(atype< Matrice_Creuse< R > * >( ), atype< pmesh3 >( ), atype< double >( ),
88                         atype< E_Array >( ));
89   }
90 
f(const basicAC_F0 & args)91   static E_F0 *f(const basicAC_F0 &args) { return new MatrixUpWind3(args); }
92 
93   AnyType operator( )(Stack s) const;
94 };
95 
gladys(double q[3][2],double u[2],double c[3],double a[3][3])96 int gladys(double q[3][2], double u[2], double c[3], double a[3][3]) {    // PSI Deconninck
97   // computes matrix a on a triangle for the Chacon-Reina Petrof-Galerkin upwind
98   // working arrays
99   double dw[3][2];                 // basis function gradients times  area
100   double ua[2], kk[3], beta[3];    // to define a[][]
101   double udc = 0;                  // u.grad(w)*area
102   bool oneaval = false;
103   int i1 = -1;
104 
105   for (int i = 0; i < 3; i++) {
106     int ip = (i + 1) % 3, ipp = (ip + 1) % 3;
107 
108     for (int j = 0; j < 2; j++) {
109       dw[i][1 - j] = (2 * j - 1) * (q[ipp][j] - q[ip][j]) / 2;
110     }
111   }
112 
113   for (int i = 0; i < 3; i++) {
114     kk[i] = u[0] * dw[i][0] + u[1] * dw[i][1];
115     udc += kk[i] * c[i];
116   }
117 
118   for (int i = 0; i < 3; i++) {
119     ua[0] = u[0];
120     ua[1] = u[1];
121     int ip = (i + 1) % 3, ipp = (ip + 1) % 3;
122     if (kk[i] > 0 && kk[ip] <= 0 && kk[ipp] <= 0) {
123       beta[i] = 1;
124       beta[ip] = 0;
125       beta[ipp] = 0;
126       oneaval = true;
127     } else if (kk[i] <= 0 && kk[ip] > 0 && kk[ipp] > 0) {
128       i1 = i;
129     }
130   }
131 
132   if (!oneaval) {
133     if (i1 < 0) {
134       cout << "bug\n";
135     }
136 
137     int i = i1, ip = (i + 1) % 3, ipp = (i + 2) % 3;
138     double lambda = (c[ip] - c[i]) * (c[ipp] - c[i]);
139     if (fabs(lambda) < -1e-20) {
140       return 0;
141     }
142 
143     if (lambda < 0) {
144       if (udc > 0) {
145         beta[i] = 0;
146         beta[ip] = 0;
147         beta[ipp] = 1;
148         ua[0] = udc * (q[ipp][0] - q[i][0]) / (c[ipp] - c[i]);
149         ua[1] = udc * (q[ipp][1] - q[i][1]) / (c[ipp] - c[i]);
150       } else {
151         beta[i] = 0;
152         beta[ipp] = 0;
153         beta[ip] = 1;
154         ua[0] = udc * (q[ip][0] - q[i][0]) / (c[ip] - c[i]);
155         ua[1] = udc * (q[ip][1] - q[i][1]) / (c[ip] - c[i]);
156       }
157     } else {
158       beta[i] = 0;
159       beta[ip] = kk[ip] * (c[ip] - c[i]) / udc;
160       beta[ipp] = kk[ipp] * (c[ipp] - c[i]) / udc;
161     }
162   }
163 
164   for (int i = 0; i < 3; i++) {
165     for (int j = 0; j < 3; j++) {
166       a[i][j] = beta[i] * (ua[0] * dw[j][0] + ua[1] * dw[j][1]);
167     }
168   }
169 
170   return 1;
171 }
172 
operator ( )(Stack stack) const173 AnyType MatrixUpWind0::operator( )(Stack stack) const {
174   Matrice_Creuse< R > *sparce_mat = GetAny< Matrice_Creuse< R > * >((*emat)(stack));
175   MatriceMorse< R > *amorse = 0;
176   MeshPoint *mp(MeshPointStack(stack)), mps = *mp;
177   const Mesh *pTh = GetAny< pmesh >((*expTh)(stack));
178   ffassert(pTh);
179   const Mesh &Th(*pTh);
180   {
181     MatriceMorse< R > *pAij = new MatriceMorse< R >(Th.nv,-1,0,0), &Aij = *pAij;
182 
183     KN< double > cc(Th.nv);
184     double infini = DBL_MAX;
185     cc = infini;
186 
187     for (int it = 0; it < Th.nt; it++) {
188       for (int iv = 0; iv < 3; iv++) {
189         int i = Th(it, iv);
190         if (cc[i] == infini) {    // if nuset the set
191           mp->setP(&Th, it, iv);
192           cc[i] = GetAny< double >((*expc)(stack));
193         }
194       }
195     }
196 
197     for (int k = 0; k < Th.nt; k++) {
198       const Triangle &K(Th[k]);
199       const Vertex &A(K[0]), &B(K[1]), &C(K[2]);
200       R2 Pt(1. / 3., 1. / 3.);
201       R u[2];
202       MeshPointStack(stack)->set(Th, K(Pt), Pt, K, K.lab);
203       u[0] = GetAny< R >((*expu1)(stack));
204       u[1] = GetAny< R >((*expu2)(stack));
205 
206       int ii[3] = {Th(A), Th(B), Th(C)};
207       double q[3][2] = {{A.x, A.y}, {B.x, B.y}, {C.x, C.y}};    // coordinates of 3 vertices (input)
208       double c[3] = {cc[ii[0]], cc[ii[1]], cc[ii[2]]};
209       double a[3][3];
210       if (gladys(q, u, c, a)) {
211         for (int i = 0; i < 3; i++) {
212           for (int j = 0; j < 3; j++) {
213             if (fabs(a[i][j]) >= 1e-30) {
214               Aij[make_pair(ii[i], ii[j])] += a[i][j];
215             }
216           }
217         }
218       }
219     }
220 
221     amorse = pAij;    // new MatriceMorse<R>(Th.nv, Th.nv, Aij, false);
222   }
223   sparce_mat->Uh = UniqueffId( );
224   sparce_mat->Vh = UniqueffId( );
225   sparce_mat->A.master(amorse);
226   sparce_mat->typemat =
227     0;    //(amorse->n == amorse->m) ? TypeSolveMat(TypeSolveMat::GMRES) :
228           // TypeSolveMat(TypeSolveMat::NONESQUARE);// none square matrice (morse)
229   *mp = mps;
230 
231   if (verbosity > 3) {
232     cout << "  End Build MatrixUpWind : " << endl;
233   }
234 
235   return sparce_mat;
236 }
237 
Marco(const Mesh3::Element & K,R3 U,R c[4],double a[4][4])238 int Marco(const Mesh3::Element &K, R3 U, R c[4], double a[4][4]) {    // PSI Deconninck
239   ExecError("Not Implemented Sorry Marco!");
240   return 0;
241 }
242 
operator ( )(Stack stack) const243 AnyType MatrixUpWind3::operator( )(Stack stack) const {
244   Matrice_Creuse< R > *sparce_mat = GetAny< Matrice_Creuse< R > * >((*emat)(stack));
245   MatriceMorse< R > *amorse = 0;
246   MeshPoint *mp(MeshPointStack(stack)), mps = *mp;
247   const Mesh3 *pTh = GetAny< pmesh3 >((*expTh)(stack));
248   ffassert(pTh);
249   const Mesh3 &Th(*pTh);
250   {
251     MatriceMorse< R > *pAij = new MatriceMorse< R >(Th.nv,-1,0,0), &Aij = *pAij;
252 
253     KN< double > cc(Th.nv);
254     double infini = DBL_MAX;
255     cc = infini;
256 
257     for (int it = 0; it < Th.nt; it++) {
258       for (int iv = 0; iv < 4; iv++) {
259         int i = Th(it, iv);
260         if (cc[i] == infini) {    // if nuset the set
261           mp->setP(&Th, it, iv);
262           cc[i] = GetAny< double >((*expc)(stack));
263         }
264       }
265     }
266 
267     for (int k = 0; k < Th.nt; k++) {
268       const Mesh3::Element &K(Th[k]);
269       const Mesh3::Vertex &A(K[0]), &B(K[1]), &C(K[2]), &D(K[3]);
270       R3 Pt(1. / 4., 1. / 4., 1. / 4.);
271       R3 U;
272       MeshPointStack(stack)->set(Th, K(Pt), Pt, K, K.lab);
273       U.x = GetAny< R >((*expu1)(stack));
274       U.y = GetAny< R >((*expu2)(stack));
275       U.z = GetAny< R >((*expu3)(stack));
276 
277       int ii[4] = {Th(A), Th(B), Th(C), Th(D)};    // number of 4 vertex
278       double c[4] = {cc[ii[0]], cc[ii[1]], cc[ii[2]], cc[ii[3]]};
279       double a[4][4];
280       if (Marco(K, U, c, a)) {
281         for (int i = 0; i < 4; i++) {
282           for (int j = 0; j < 4; j++) {
283             if (fabs(a[i][j]) >= 1e-30) {
284               Aij[make_pair(ii[i], ii[j])] += a[i][j];
285             }
286           }
287         }
288       }
289     }
290 
291     amorse = pAij;    // new MatriceMorse<R>(Th.nv, Th.nv, Aij, false);
292   }
293   sparce_mat->Uh = UniqueffId( );
294   sparce_mat->Vh = UniqueffId( );
295   sparce_mat->A.master(amorse);
296   sparce_mat->typemat =
297     0;    //(amorse->n == amorse->m) ? TypeSolveMat(TypeSolveMat::GMRES) :
298           // TypeSolveMat(TypeSolveMat::NONESQUARE);// none square matrice (morse)
299   *mp = mps;
300 
301   if (verbosity > 3) {
302     cout << "  End Build MatrixUpWind : " << endl;
303   }
304 
305   return sparce_mat;
306 }
307 
Load_Init()308 static void Load_Init( ) {
309   cout << " lood: init Mat Chacon " << endl;
310   Global.Add("MatUpWind0", "(", new OneOperatorCode< MatrixUpWind0 >( ));
311   Global.Add("MatUpWind0", "(", new OneOperatorCode< MatrixUpWind3 >( ));
312 }
313 
314 LOADFUNC(Load_Init)
315