1 // GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege
2 //
3 // See the LICENSE.txt file for license information. Please report all
4 // issues on https://gitlab.onelab.info/getdp/getdp/issues.
5
6 #include "GetDPConfig.h"
7 #include "ProData.h"
8 #include "F.h"
9 #include "Message.h"
10
11 extern struct CurrentData Current ;
12
13 #if defined(HAVE_GMSH)
14
15 #include <gmsh/GmshGlobal.h>
16 #include <gmsh/PView.h>
17 #include <gmsh/PViewData.h>
18 #include <gmsh/Numeric.h>
19
F_Field(F_ARG)20 void F_Field(F_ARG)
21 {
22 if(A->Type != VECTOR){
23 Message::Error("Field[] expects XYZ coordinates as argument");
24 return;
25 }
26 if(PView::list.empty()){
27 Message::Error("No views available to interpolate from");
28 return;
29 }
30 double x = A->Val[0];
31 double y = A->Val[1];
32 double z = A->Val[2];
33
34 if(Fct->NbrArguments > 1){
35 Message::Error("Time and additional arguments are not supported in Field: "
36 "use {Scalar,Vector,Tensor}Field instead");
37 return;
38 }
39
40 for (int k = 0; k < Current.NbrHar; k++)
41 for (int j = 0; j < 9; j++)
42 V->Val[MAX_DIM * k + j] = 0. ;
43 V->Type = SCALAR;
44
45 std::vector<int> iview;
46 if(!Fct->NbrParameters){
47 // use last view by default
48 iview.push_back(PView::list.back()->getTag());
49 }
50 else{
51 for(int i = 0; i < Fct->NbrParameters; i++)
52 iview.push_back(Fct->Para[i]);
53 }
54
55 double N = 0.;
56
57 // add the values from all specified views
58 for(unsigned int i = 0; i < iview.size(); i++){
59
60 PView *v = PView::getViewByTag(iview[i]);
61 if(!v){
62 Message::Error("View with tag %d does not exist", iview[i]);
63 return;
64 }
65
66 PViewData *data = v->getData();
67
68 std::vector<double> val(9 * data->getNumTimeSteps());
69 if(data->searchScalar(x, y, z, &val[0])){
70 V->Val[0] += val[0];
71 if(Current.NbrHar == 2 && data->getNumTimeSteps() > 1)
72 V->Val[MAX_DIM] += val[1];
73 V->Type = SCALAR;
74 N += 1.;
75 }
76 else if(data->searchVector(x, y, z, &val[0])){
77 for(int j = 0; j < 3; j++)
78 V->Val[j] += val[j];
79 if(Current.NbrHar == 2 && data->getNumTimeSteps() > 1){
80 for(int j = 0; j < 3; j++)
81 V->Val[MAX_DIM + j] += val[3 + j];
82 }
83 V->Type = VECTOR;
84 N += 1.;
85 }
86 else if(data->searchTensor(x, y, z, &val[0])){
87 for(int j = 0; j < 9; j++)
88 V->Val[j] += val[j];
89 if(Current.NbrHar == 2 && data->getNumTimeSteps() > 1){
90 for(int j = 0; j < 9; j++)
91 V->Val[MAX_DIM + j] += val[9 + j];
92 }
93 V->Type = TENSOR;
94 N += 1.;
95 }
96 else{
97 Message::Error("Did not find data at point (%g,%g,%g) in View with tag %d",
98 x, y, z, iview[i]);
99 }
100 }
101
102 if(N > 1.){
103 Message::Debug("Averaging data %g times on vertex (%g,%g,%g)", N, x, y, z);
104 for (int k = 0; k < Current.NbrHar; k++)
105 for (int j = 0; j < 9; j++)
106 V->Val[MAX_DIM * k + j] = 0. ;
107 }
108 }
109
F_X_Field(F_ARG,int type,bool complex,bool grad=false)110 static void F_X_Field(F_ARG, int type, bool complex, bool grad=false)
111 {
112 if(A->Type != VECTOR){
113 Message::Error("Field[] expects XYZ coordinates as argument");
114 return;
115 }
116 if(PView::list.empty()){
117 Message::Error("No views available to interpolate from");
118 return;
119 }
120
121 double x = A->Val[0];
122 double y = A->Val[1];
123 double z = A->Val[2];
124 int numComp = (type == SCALAR) ? (grad ? 3 : 1) :
125 (type == VECTOR) ? (grad ? 9 : 3) : 9; // TODO: grad of tensor
126
127 int NbrArg = Fct->NbrArguments ;
128 int TimeStep = 0, MatchElement = 0;
129 if(NbrArg >= 2){
130 if((A+1)->Type != SCALAR){
131 Message::Error("Expected scalar second argument (time step)");
132 return;
133 }
134 TimeStep = (int)(A+1)->Val[0];
135 }
136 if(NbrArg >= 3){
137 if((A+2)->Type != SCALAR){
138 Message::Error("Expected scalar second argument (element matching flag)");
139 return;
140 }
141 MatchElement = (int)(A+2)->Val[0];
142 }
143
144 // TODO: we could treat the third arguement as a tolerance (and call
145 // searchScalarWithTol & friends)
146
147 // Complex{Scalar,Vector,Tensor}Field assume that the Gmsh view contains real
148 // and imaginary parts for each step
149 if(complex) TimeStep *= 2;
150
151 for (int k = 0; k < Current.NbrHar; k++)
152 for (int j = 0; j < numComp; j++)
153 V->Val[MAX_DIM * k + j] = 0. ;
154 V->Type = (numComp == 1) ? SCALAR : (numComp == 3) ? VECTOR : TENSOR;
155
156 std::vector<int> iview;
157 if(!Fct->NbrParameters){
158 // use last view by default
159 iview.push_back(PView::list.back()->getTag());
160 }
161 else{
162 for(int i = 0; i < Fct->NbrParameters; i++)
163 iview.push_back(Fct->Para[i]);
164 }
165
166 int qn = 0;
167 double *qx = 0, *qy = 0, *qz = 0;
168 if(Current.Element){
169 qn = MatchElement ? Current.Element->GeoElement->NbrNodes : 0;
170 qx = Current.Element->x;
171 qy = Current.Element->y;
172 qz = Current.Element->z;
173 }
174
175 double N = 0.;
176
177 // add the values from all specified views
178 for(unsigned int i = 0; i < iview.size(); i++){
179 PView *v = PView::getViewByTag(iview[i]);
180 if(!v){
181 Message::Error("View with tag %d does not exist", iview[i]);
182 return;
183 }
184 PViewData *data = v->getData();
185 if(TimeStep < 0 || TimeStep >= data->getNumTimeSteps()){
186 Message::Error("Invalid step %d in View with tag %d", TimeStep, iview[i]);
187 continue;
188 }
189 std::vector<double> val(numComp * data->getNumTimeSteps());
190 bool found = false;
191 switch(type){
192 case SCALAR :
193 if(data->searchScalar(x, y, z, &val[0], -1, 0, qn, qx, qy, qz, grad))
194 found = true;
195 break;
196 case VECTOR :
197 if(data->searchVector(x, y, z, &val[0], -1, 0, qn, qx, qy, qz, grad))
198 found = true;
199 break;
200 case TENSOR :
201 // TODO: grad of tensor not allowed yet - not sure how we should return
202 // the values; provide 3 routines that return 3 tensors, or add argumemt
203 // to select what to return?
204 if(data->searchTensor(x, y, z, &val[0], -1, 0, qn, qx, qy, qz, false))
205 found = true;
206 break;
207 }
208 if(found){
209 for(int j = 0; j < numComp; j++)
210 V->Val[j] += val[numComp * TimeStep + j];
211 if(complex && Current.NbrHar == 2 && data->getNumTimeSteps() > TimeStep + 1)
212 for(int j = 0; j < numComp; j++)
213 V->Val[MAX_DIM + j] += val[numComp * (TimeStep + 1) + j];
214 N += 1.;
215 }
216 }
217
218 if(N > 1.){
219 Message::Debug("Averaging data %g times on vertex (%g,%g,%g)", N, x, y, z);
220 for (int k = 0; k < Current.NbrHar; k++)
221 for (int j = 0; j < numComp; j++)
222 V->Val[MAX_DIM * k + j] /= N ;
223 }
224 }
225
F_Distance(F_ARG)226 void F_Distance(F_ARG)
227 {
228 if(A->Type != VECTOR){
229 Message::Error("Distance[] expects XYZ coordinates as argument");
230 return;
231 }
232
233 SPoint3 p(A->Val[0], A->Val[1], A->Val[2]);
234
235 if(Fct->NbrParameters != 1){
236 Message::Error("Distance[] requires one view tag as parameter");
237 return;
238 }
239
240 PView *v = PView::getViewByTag(Fct->Para[0]);
241 if(!v){
242 Message::Error("View with tag %d does not exist", Fct->Para[0]);
243 return;
244 }
245
246 PViewData *data = v->getData();
247 if(data->getNumTimeSteps() != 1 || !data->hasTimeStep(0)) {
248 Message::Error("View used in Distance[] should have a single step");
249 return;
250 }
251
252 double dist = 1e200;
253
254 for(int ent = 0; ent < data->getNumEntities(0); ent++) {
255 if(data->skipEntity(0, ent)) continue;
256 for(int ele = 0; ele < data->getNumElements(0, ent); ele++) {
257 if(data->skipElement(0, ent, ele)) continue;
258 int numComp = data->getNumComponents(0, ent, ele);
259 if(numComp != 3) {
260 Message::Error("Distance[] expects a vector-valued view");
261 return;
262 }
263 int numNodes = data->getNumNodes(0, ent, ele);
264 if(numNodes != 2 && numNodes != 3) {
265 Message::Error("Distance[] can only currently handle lines and triangles");
266 return;
267 }
268 double nx[3], ny[3], nz[3], dx[3], dy[3], dz[3];
269 for(int nod = 0; nod < numNodes; nod++) {
270 data->getValue(0, ent, ele, nod, 0, dx[nod]);
271 data->getValue(0, ent, ele, nod, 1, dy[nod]);
272 data->getValue(0, ent, ele, nod, 2, dz[nod]);
273 data->getNode(0, ent, ele, nod, nx[nod], ny[nod], nz[nod]);
274 }
275 double d;
276 SPoint3 p1(nx[0] + dx[0], ny[0] + dy[0], nz[0] + dz[0]);
277 SPoint3 p2(nx[1] + dx[1], ny[1] + dy[1], nz[1] + dz[1]);
278 SPoint3 closePt;
279 if(numNodes == 2) {
280 // Warning: this is currently NOT signed!
281 signedDistancePointLine(p1, p2, p, d, closePt);
282 }
283 else if(numNodes == 3) {
284 SPoint3 p3(nx[2] + dx[2], ny[2] + dy[2], nz[2] + dz[2]);
285 signedDistancePointTriangle(p1, p2, p3, p, d, closePt);
286 }
287 dist = std::min(dist, d);
288 }
289 }
290
291 V->Type = SCALAR ;
292 V->Val[0] = dist;
293 V->Val[MAX_DIM] = 0.;
294 for (int k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) {
295 V->Val[MAX_DIM* k] = V->Val[0] ;
296 V->Val[MAX_DIM* (k+1)] = 0. ;
297 }
298 }
299
300 #else
301
F_Field(F_ARG)302 void F_Field(F_ARG)
303 {
304 Message::Error("You need to compile GetDP with Gmsh support to use Field[]");
305 V->Val[0] = 0. ;
306 V->Type = SCALAR ;
307 }
308
F_X_Field(F_ARG,int type,bool complex,bool grad=false)309 static void F_X_Field(F_ARG, int type, bool complex, bool grad=false)
310 {
311 Message::Error("You need to compile GetDP with Gmsh support to use Field[]");
312 V->Val[0] = 0. ;
313 V->Type = SCALAR ;
314 }
315
F_Distance(F_ARG)316 void F_Distance(F_ARG)
317 {
318 Message::Error("You need to compile GetDP with Gmsh support to use Distance[]");
319 V->Val[0] = 0. ;
320 V->Type = SCALAR ;
321 }
322
323 #endif
324
F_ScalarField(F_ARG)325 void F_ScalarField(F_ARG){ F_X_Field(Fct, A, V, SCALAR, false); }
F_VectorField(F_ARG)326 void F_VectorField(F_ARG){ F_X_Field(Fct, A, V, VECTOR, false); }
F_TensorField(F_ARG)327 void F_TensorField(F_ARG){ F_X_Field(Fct, A, V, TENSOR, false); }
F_ComplexScalarField(F_ARG)328 void F_ComplexScalarField(F_ARG){ F_X_Field(Fct, A, V, SCALAR, true); }
F_ComplexVectorField(F_ARG)329 void F_ComplexVectorField(F_ARG){ F_X_Field(Fct, A, V, VECTOR, true); }
F_ComplexTensorField(F_ARG)330 void F_ComplexTensorField(F_ARG){ F_X_Field(Fct, A, V, TENSOR, true); }
331
F_GradScalarField(F_ARG)332 void F_GradScalarField(F_ARG){ F_X_Field(Fct, A, V, SCALAR, false, true); }
F_GradVectorField(F_ARG)333 void F_GradVectorField(F_ARG){ F_X_Field(Fct, A, V, VECTOR, false, true); }
F_GradComplexScalarField(F_ARG)334 void F_GradComplexScalarField(F_ARG){ F_X_Field(Fct, A, V, SCALAR, true, true); }
F_GradComplexVectorField(F_ARG)335 void F_GradComplexVectorField(F_ARG){ F_X_Field(Fct, A, V, VECTOR, true, true); }
336