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