Nektar++
ElUtil.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ElUtil.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Calculate jacobians of elements.
32 //
33 ////////////////////////////////////////////////////////////////////////////////
34 
35 #include "ElUtil.h"
36 #include "ProcessVarOpti.h"
37 
38 #include <mutex>
39 #include <LibUtilities/Foundations/ManagerAccess.h>
40 
41 using namespace std;
42 
43 namespace Nektar
44 {
45 namespace Utilities
46 {
47 
48 std::mutex mtx2;
49 
50 ElUtil::ElUtil(ElementSharedPtr e, DerivUtilSharedPtr d, ResidualSharedPtr r,
51  int n, int o)
52 {
53  m_el = e;
54  m_derivUtil = d;
55  m_res = r;
56  m_mode = n;
57  m_order = o;
58  m_dim = m_el->GetDim();
59  vector<NodeSharedPtr> ns;
60  m_el->GetCurvedNodes(ns);
61  nodes.resize(ns.size());
62  for (int i = 0; i < ns.size(); ++i)
63  {
64  nodes[i].resize(m_dim);
65  nodes[i][0] = &ns[i]->m_x;
66 
67  if (m_dim >= 2)
68  {
69  nodes[i][1] = &ns[i]->m_y;
70  }
71 
72  if (m_dim >= 3)
73  {
74  nodes[i][2] = &ns[i]->m_z;
75  }
76 
77  m_idmap[ns[i]->m_id] = i;
78  }
79  MappingIdealToRef();
80 }
81 
82 void ElUtil::MappingIdealToRef()
83 {
84  if (m_el->GetConf().m_e == LibUtilities::eQuadrilateral)
85  {
86  LibUtilities::PointsKey pkey1(m_mode, LibUtilities::eNodalQuadElec);
87  LibUtilities::PointsKey pkey2(m_mode + m_order,
88  LibUtilities::eNodalQuadElec);
89 
90  Array<OneD, NekDouble> u1(m_derivUtil->ptsStd), v1(m_derivUtil->ptsStd),
91  u2(m_derivUtil->pts), v2(m_derivUtil->pts);
92 
93  LibUtilities::PointsManager()[pkey1]->GetPoints(u1, v1);
94  LibUtilities::PointsManager()[pkey2]->GetPoints(u2, v2);
95 
96  vector<vector<NekDouble> > xyz(4);
97  vector<NodeSharedPtr> ns = m_el->GetVertexList();
98  for (int i = 0; i < 4; i++)
99  {
100  vector<NekDouble> x(3);
101  x[0] = ns[i]->m_x;
102  x[1] = ns[i]->m_y;
103  x[2] = ns[i]->m_z;
104  xyz[i] = x;
105  }
106 
107  for (int i = 0; i < m_derivUtil->ptsStd; ++i)
108  {
109  NekDouble a1 = 0.5 * (1 - u1[i]);
110  NekDouble a2 = 0.5 * (1 + u1[i]);
111  NekDouble b1 = 0.5 * (1 - v1[i]);
112  NekDouble b2 = 0.5 * (1 + v1[i]);
113 
114  DNekMat J(2, 2, 1.0, eFULL);
115 
116  J(0, 0) = -0.5 * b1 * xyz[0][0] + 0.5 * b1 * xyz[1][0] +
117  0.5 * b2 * xyz[2][0] - 0.5 * b2 * xyz[3][0];
118  J(1, 0) = -0.5 * b1 * xyz[0][1] + 0.5 * b1 * xyz[1][1] +
119  0.5 * b2 * xyz[2][1] - 0.5 * b2 * xyz[3][1];
120 
121  J(0, 1) = -0.5 * a1 * xyz[0][0] - 0.5 * a2 * xyz[1][0] +
122  0.5 * a2 * xyz[2][0] + 0.5 * a1 * xyz[3][0];
123  J(1, 1) = -0.5 * a1 * xyz[0][1] - 0.5 * a2 * xyz[1][1] +
124  0.5 * a2 * xyz[2][1] + 0.5 * a1 * xyz[3][1];
125 
126  J.Invert();
127 
128  vector<NekDouble> r(10, 0.0); // store det in 10th entry
129 
130  r[9] = 1.0 / (J(0, 0) * J(1, 1) - J(0, 1) * J(1, 0));
131 
132  r[0] = J(0, 0);
133  r[1] = J(1, 0);
134  r[2] = 0.0;
135  r[3] = J(0, 1);
136  r[4] = J(1, 1);
137  r[5] = 0.0;
138  r[6] = 0.0;
139  r[7] = 0.0;
140  r[8] = 0.0;
141  mapsStd.push_back(r);
142  }
143 
144  for (int i = 0; i < m_derivUtil->pts; ++i)
145  {
146  NekDouble a1 = 0.5 * (1 - u2[i]);
147  NekDouble a2 = 0.5 * (1 + u2[i]);
148  NekDouble b1 = 0.5 * (1 - v2[i]);
149  NekDouble b2 = 0.5 * (1 + v2[i]);
150 
151  DNekMat J(2, 2, 1.0, eFULL);
152 
153  J(0, 0) = -0.5 * b1 * xyz[0][0] + 0.5 * b1 * xyz[1][0] +
154  0.5 * b2 * xyz[2][0] - 0.5 * b2 * xyz[3][0];
155  J(1, 0) = -0.5 * b1 * xyz[0][1] + 0.5 * b1 * xyz[1][1] +
156  0.5 * b2 * xyz[2][1] - 0.5 * b2 * xyz[3][1];
157 
158  J(0, 1) = -0.5 * a1 * xyz[0][0] - 0.5 * a2 * xyz[1][0] +
159  0.5 * a2 * xyz[2][0] + 0.5 * a1 * xyz[3][0];
160  J(1, 1) = -0.5 * a1 * xyz[0][1] - 0.5 * a2 * xyz[1][1] +
161  0.5 * a2 * xyz[2][1] + 0.5 * a1 * xyz[3][1];
162 
163  J.Invert();
164 
165  vector<NekDouble> r(10, 0.0); // store det in 10th entry
166 
167  r[9] = 1.0 / (J(0, 0) * J(1, 1) - J(0, 1) * J(1, 0));
168 
169  r[0] = J(0, 0);
170  r[1] = J(1, 0);
171  r[2] = 0.0;
172  r[3] = J(0, 1);
173  r[4] = J(1, 1);
174  r[5] = 0.0;
175  r[6] = 0.0;
176  r[7] = 0.0;
177  r[8] = 0.0;
178  maps.push_back(r);
179  }
180  }
181  else if (m_el->GetConf().m_e == LibUtilities::eTriangle)
182  {
183  DNekMat J(2, 2, 0.0);
184  J(0, 0) = (*nodes[1][0] - *nodes[0][0]);
185  J(1, 0) = (*nodes[1][1] - *nodes[0][1]);
186  J(0, 1) = (*nodes[2][0] - *nodes[0][0]);
187  J(1, 1) = (*nodes[2][1] - *nodes[0][1]);
188 
189  J.Invert();
190 
191  DNekMat R(2, 2, 0.0);
192  R(0, 0) = 2.0;
193  R(1, 1) = 2.0;
194 
195  J = J * R;
196 
197  for (int i = 0; i < m_derivUtil->pts; i++)
198  {
199  vector<NekDouble> r(10, 0.0); // store det in 10th entry
200 
201  r[9] = 1.0 / (J(0, 0) * J(1, 1) - J(0, 1) * J(1, 0));
202  r[0] = J(0, 0);
203  r[1] = J(1, 0);
204  r[2] = 0.0;
205  r[3] = J(0, 1);
206  r[4] = J(1, 1);
207  r[5] = 0.0;
208  r[6] = 0.0;
209  r[7] = 0.0;
210  r[8] = 0.0;
211  maps.push_back(r);
212  mapsStd.push_back(r);
213  }
214  }
215  else if (m_el->GetConf().m_e == LibUtilities::eTetrahedron)
216  {
217  DNekMat J(3, 3, 0.0);
218  J(0, 0) = (*nodes[1][0] - *nodes[0][0]);
219  J(1, 0) = (*nodes[1][1] - *nodes[0][1]);
220  J(2, 0) = (*nodes[1][2] - *nodes[0][2]);
221  J(0, 1) = (*nodes[2][0] - *nodes[0][0]);
222  J(1, 1) = (*nodes[2][1] - *nodes[0][1]);
223  J(2, 1) = (*nodes[2][2] - *nodes[0][2]);
224  J(0, 2) = (*nodes[3][0] - *nodes[0][0]);
225  J(1, 2) = (*nodes[3][1] - *nodes[0][1]);
226  J(2, 2) = (*nodes[3][2] - *nodes[0][2]);
227 
228  J.Invert();
229 
230  DNekMat R(3, 3, 0.0);
231  R(0, 0) = 2.0;
232  R(1, 1) = 2.0;
233  R(2, 2) = 2.0;
234 
235  J = J * R;
236 
237  for (int i = 0; i < m_derivUtil->pts; i++)
238  {
239  vector<NekDouble> r(10, 0.0); // store det in 10th entry
240 
241  r[9] = 1.0 / (J(0, 0) * (J(1, 1) * J(2, 2) - J(2, 1) * J(1, 2)) -
242  J(0, 1) * (J(1, 0) * J(2, 2) - J(2, 0) * J(1, 2)) +
243  J(0, 2) * (J(1, 0) * J(2, 1) - J(2, 0) * J(1, 1)));
244 
245  r[0] = J(0, 0);
246  r[1] = J(1, 0);
247  r[2] = J(2, 0);
248  r[3] = J(0, 1);
249  r[4] = J(1, 1);
250  r[5] = J(2, 1);
251  r[6] = J(0, 2);
252  r[7] = J(1, 2);
253  r[8] = J(2, 2);
254  maps.push_back(r);
255  mapsStd.push_back(r);
256  }
257  }
258  else if (m_el->GetConf().m_e == LibUtilities::ePrism)
259  {
260  LibUtilities::PointsKey pkey1(m_mode, LibUtilities::eNodalPrismElec);
261  LibUtilities::PointsKey pkey2(m_mode + m_order,
262  LibUtilities::eNodalPrismSPI);
263  Array<OneD, NekDouble> u1, v1, u2, v2, w1, w2;
264  LibUtilities::PointsManager()[pkey1]->GetPoints(u1, v1, w1);
265  LibUtilities::PointsManager()[pkey2]->GetPoints(u2, v2, w2);
266 
267  vector<vector<NekDouble> > xyz(6);
268  vector<NodeSharedPtr> ns = m_el->GetVertexList();
269  for (int i = 0; i < 6; i++)
270  {
271  vector<NekDouble> x(3);
272  x[0] = ns[i]->m_x;
273  x[1] = ns[i]->m_y;
274  x[2] = ns[i]->m_z;
275  xyz[i] = x;
276  }
277 
278  for (int i = 0; i < m_derivUtil->ptsStd; ++i)
279  {
280  NekDouble a2 = 0.5 * (1 + u1[i]);
281  NekDouble b1 = 0.5 * (1 - v1[i]);
282  NekDouble b2 = 0.5 * (1 + v1[i]);
283  NekDouble c2 = 0.5 * (1 + w1[i]);
284  NekDouble d = 0.5 * (u1[i] + w1[i]);
285 
286  DNekMat J(3, 3, 1.0, eFULL);
287 
288  J(0, 0) = -0.5 * b1 * xyz[0][0] + 0.5 * b1 * xyz[1][0] +
289  0.5 * b2 * xyz[2][0] - 0.5 * b2 * xyz[3][0];
290  J(1, 0) = -0.5 * b1 * xyz[0][1] + 0.5 * b1 * xyz[1][1] +
291  0.5 * b2 * xyz[2][1] - 0.5 * b2 * xyz[3][1];
292  J(2, 0) = -0.5 * b1 * xyz[0][2] + 0.5 * b1 * xyz[1][2] +
293  0.5 * b2 * xyz[2][2] - 0.5 * b2 * xyz[3][2];
294 
295  J(0, 1) = 0.5 * d * xyz[0][0] - 0.5 * a2 * xyz[1][0] +
296  0.5 * a2 * xyz[2][0] - 0.5 * d * xyz[3][0] -
297  0.5 * c2 * xyz[4][0] + 0.5 * c2 * xyz[5][0];
298  J(1, 1) = 0.5 * d * xyz[0][1] - 0.5 * a2 * xyz[1][1] +
299  0.5 * a2 * xyz[2][1] - 0.5 * d * xyz[3][1] -
300  0.5 * c2 * xyz[4][1] + 0.5 * c2 * xyz[5][1];
301  J(2, 1) = 0.5 * d * xyz[0][2] - 0.5 * a2 * xyz[1][2] +
302  0.5 * a2 * xyz[2][2] - 0.5 * d * xyz[3][2] -
303  0.5 * c2 * xyz[4][2] + 0.5 * c2 * xyz[5][2];
304 
305  J(0, 2) = -0.5 * b1 * xyz[0][0] - 0.5 * b2 * xyz[3][0] +
306  0.5 * b1 * xyz[4][0] + 0.5 * b2 * xyz[5][0];
307  J(1, 2) = -0.5 * b1 * xyz[0][1] - 0.5 * b2 * xyz[3][1] +
308  0.5 * b1 * xyz[4][1] + 0.5 * b2 * xyz[5][1];
309  J(2, 2) = -0.5 * b1 * xyz[0][2] - 0.5 * b2 * xyz[3][2] +
310  0.5 * b1 * xyz[4][2] + 0.5 * b2 * xyz[5][2];
311 
312  J.Invert();
313 
314  vector<NekDouble> r(10, 0.0); // store det in 10th entry
315 
316  r[9] = 1.0 / (J(0, 0) * (J(1, 1) * J(2, 2) - J(2, 1) * J(1, 2)) -
317  J(0, 1) * (J(1, 0) * J(2, 2) - J(2, 0) * J(1, 2)) +
318  J(0, 2) * (J(1, 0) * J(2, 1) - J(2, 0) * J(1, 1)));
319 
320  r[0] = J(0, 0);
321  r[1] = J(1, 0);
322  r[2] = J(2, 0);
323  r[3] = J(0, 1);
324  r[4] = J(1, 1);
325  r[5] = J(2, 1);
326  r[6] = J(0, 2);
327  r[7] = J(1, 2);
328  r[8] = J(2, 2);
329  mapsStd.push_back(r);
330  }
331  for (int i = 0; i < m_derivUtil->pts; ++i)
332  {
333  NekDouble a2 = 0.5 * (1 + u2[i]);
334  NekDouble b1 = 0.5 * (1 - v2[i]);
335  NekDouble b2 = 0.5 * (1 + v2[i]);
336  NekDouble c2 = 0.5 * (1 + w2[i]);
337  NekDouble d = 0.5 * (u2[i] + w2[i]);
338 
339  DNekMat J(3, 3, 1.0, eFULL);
340 
341  J(0, 0) = -0.5 * b1 * xyz[0][0] + 0.5 * b1 * xyz[1][0] +
342  0.5 * b2 * xyz[2][0] - 0.5 * b2 * xyz[3][0];
343  J(1, 0) = -0.5 * b1 * xyz[0][1] + 0.5 * b1 * xyz[1][1] +
344  0.5 * b2 * xyz[2][1] - 0.5 * b2 * xyz[3][1];
345  J(2, 0) = -0.5 * b1 * xyz[0][2] + 0.5 * b1 * xyz[1][2] +
346  0.5 * b2 * xyz[2][2] - 0.5 * b2 * xyz[3][2];
347 
348  J(0, 1) = 0.5 * d * xyz[0][0] - 0.5 * a2 * xyz[1][0] +
349  0.5 * a2 * xyz[2][0] - 0.5 * d * xyz[3][0] -
350  0.5 * c2 * xyz[4][0] + 0.5 * c2 * xyz[5][0];
351  J(1, 1) = 0.5 * d * xyz[0][1] - 0.5 * a2 * xyz[1][1] +
352  0.5 * a2 * xyz[2][1] - 0.5 * d * xyz[3][1] -
353  0.5 * c2 * xyz[4][1] + 0.5 * c2 * xyz[5][1];
354  J(2, 1) = 0.5 * d * xyz[0][2] - 0.5 * a2 * xyz[1][2] +
355  0.5 * a2 * xyz[2][2] - 0.5 * d * xyz[3][2] -
356  0.5 * c2 * xyz[4][2] + 0.5 * c2 * xyz[5][2];
357 
358  J(0, 2) = -0.5 * b1 * xyz[0][0] - 0.5 * b2 * xyz[3][0] +
359  0.5 * b1 * xyz[4][0] + 0.5 * b2 * xyz[5][0];
360  J(1, 2) = -0.5 * b1 * xyz[0][1] - 0.5 * b2 * xyz[3][1] +
361  0.5 * b1 * xyz[4][1] + 0.5 * b2 * xyz[5][1];
362  J(2, 2) = -0.5 * b1 * xyz[0][2] - 0.5 * b2 * xyz[3][2] +
363  0.5 * b1 * xyz[4][2] + 0.5 * b2 * xyz[5][2];
364 
365  J.Invert();
366 
367  vector<NekDouble> r(10, 0.0); // store det in 10th entry
368 
369  r[9] = 1.0 / (J(0, 0) * (J(1, 1) * J(2, 2) - J(2, 1) * J(1, 2)) -
370  J(0, 1) * (J(1, 0) * J(2, 2) - J(2, 0) * J(1, 2)) +
371  J(0, 2) * (J(1, 0) * J(2, 1) - J(2, 0) * J(1, 1)));
372 
373  r[0] = J(0, 0);
374  r[1] = J(1, 0);
375  r[2] = J(2, 0);
376  r[3] = J(0, 1);
377  r[4] = J(1, 1);
378  r[5] = J(2, 1);
379  r[6] = J(0, 2);
380  r[7] = J(1, 2);
381  r[8] = J(2, 2);
382  maps.push_back(r);
383  }
384  }
385  else if (m_el->GetConf().m_e == LibUtilities::eHexahedron)
386  {
387  LibUtilities::PointsKey pkey1(m_mode, LibUtilities::eNodalHexElec);
388  LibUtilities::PointsKey pkey2(m_mode + m_order,
389  LibUtilities::eNodalHexElec);
390  Array<OneD, NekDouble> u1(m_derivUtil->ptsStd), v1(m_derivUtil->ptsStd),
391  w1(m_derivUtil->ptsStd), u2(m_derivUtil->pts), v2(m_derivUtil->pts),
392  w2(m_derivUtil->pts);
393 
394  LibUtilities::PointsManager()[pkey1]->GetPoints(u1, v1, w1);
395  LibUtilities::PointsManager()[pkey2]->GetPoints(u2, v2, w2);
396 
397  vector<vector<NekDouble> > xyz(8);
398  vector<NodeSharedPtr> ns = m_el->GetVertexList();
399  for (int i = 0; i < 8; i++)
400  {
401  vector<NekDouble> x(3);
402  x[0] = ns[i]->m_x;
403  x[1] = ns[i]->m_y;
404  x[2] = ns[i]->m_z;
405  xyz[i] = x;
406  }
407 
408  for (int i = 0; i < m_derivUtil->ptsStd; ++i)
409  {
410  NekDouble a1 = 0.5 * (1 - u1[i]);
411  NekDouble a2 = 0.5 * (1 + u1[i]);
412  NekDouble b1 = 0.5 * (1 - v1[i]);
413  NekDouble b2 = 0.5 * (1 + v1[i]);
414  NekDouble c1 = 0.5 * (1 - w1[i]);
415  NekDouble c2 = 0.5 * (1 + w1[i]);
416 
417  DNekMat J(3, 3, 1.0, eFULL);
418 
419  J(0, 0) = -0.5 * b1 * c1 * xyz[0][0] + 0.5 * b1 * c1 * xyz[1][0] +
420  0.5 * b2 * c1 * xyz[2][0] - 0.5 * b2 * c1 * xyz[3][0] -
421  0.5 * b1 * c2 * xyz[5][0] + 0.5 * b1 * c2 * xyz[5][0] +
422  0.5 * b2 * c2 * xyz[6][0] - 0.5 * b2 * c2 * xyz[7][0];
423  J(1, 0) = -0.5 * b1 * c1 * xyz[0][1] + 0.5 * b1 * c1 * xyz[1][1] +
424  0.5 * b2 * c1 * xyz[2][1] - 0.5 * b2 * c1 * xyz[3][1] -
425  0.5 * b1 * c2 * xyz[5][1] + 0.5 * b1 * c2 * xyz[5][1] +
426  0.5 * b2 * c2 * xyz[6][1] - 0.5 * b2 * c2 * xyz[7][1];
427  J(2, 0) = -0.5 * b1 * c1 * xyz[0][2] + 0.5 * b1 * c1 * xyz[1][2] +
428  0.5 * b2 * c1 * xyz[2][2] - 0.5 * b2 * c1 * xyz[3][2] -
429  0.5 * b1 * c2 * xyz[5][2] + 0.5 * b1 * c2 * xyz[5][2] +
430  0.5 * b2 * c2 * xyz[6][2] - 0.5 * b2 * c2 * xyz[7][2];
431 
432  J(0, 1) = -0.5 * a1 * c1 * xyz[0][0] - 0.5 * a2 * c1 * xyz[1][0] +
433  0.5 * a2 * c1 * xyz[2][0] + 0.5 * a1 * c1 * xyz[3][0] -
434  0.5 * a1 * c2 * xyz[5][0] - 0.5 * a2 * c2 * xyz[5][0] +
435  0.5 * a2 * c2 * xyz[6][0] + 0.5 * a1 * c2 * xyz[7][0];
436  J(1, 1) = -0.5 * a1 * c1 * xyz[0][1] - 0.5 * a2 * c1 * xyz[1][1] +
437  0.5 * a2 * c1 * xyz[2][1] + 0.5 * a1 * c1 * xyz[3][1] -
438  0.5 * a1 * c2 * xyz[5][1] - 0.5 * a2 * c2 * xyz[5][1] +
439  0.5 * a2 * c2 * xyz[6][1] + 0.5 * a1 * c2 * xyz[7][1];
440  J(2, 1) = -0.5 * a1 * c1 * xyz[0][2] - 0.5 * a2 * c1 * xyz[1][2] +
441  0.5 * a2 * c1 * xyz[2][2] + 0.5 * a1 * c1 * xyz[3][2] -
442  0.5 * a1 * c2 * xyz[5][2] - 0.5 * a2 * c2 * xyz[5][2] +
443  0.5 * a2 * c2 * xyz[6][2] + 0.5 * a1 * c2 * xyz[7][2];
444 
445  J(0, 0) = -0.5 * b1 * a1 * xyz[0][0] - 0.5 * b1 * a2 * xyz[1][0] -
446  0.5 * b2 * a2 * xyz[2][0] - 0.5 * b2 * a1 * xyz[3][0] +
447  0.5 * b1 * a1 * xyz[5][0] + 0.5 * b1 * a2 * xyz[5][0] +
448  0.5 * b2 * a2 * xyz[6][0] + 0.5 * b2 * a1 * xyz[7][0];
449  J(1, 0) = -0.5 * b1 * a1 * xyz[0][1] - 0.5 * b1 * a2 * xyz[1][1] -
450  0.5 * b2 * a2 * xyz[2][1] - 0.5 * b2 * a1 * xyz[3][1] +
451  0.5 * b1 * a1 * xyz[5][1] + 0.5 * b1 * a2 * xyz[5][1] +
452  0.5 * b2 * a2 * xyz[6][1] + 0.5 * b2 * a1 * xyz[7][1];
453  J(2, 0) = -0.5 * b1 * a1 * xyz[0][2] - 0.5 * b1 * a2 * xyz[1][2] -
454  0.5 * b2 * a2 * xyz[2][2] - 0.5 * b2 * a1 * xyz[3][2] +
455  0.5 * b1 * a1 * xyz[5][2] + 0.5 * b1 * a2 * xyz[5][2] +
456  0.5 * b2 * a2 * xyz[6][2] + 0.5 * b2 * a1 * xyz[7][2];
457 
458  J.Invert();
459 
460  vector<NekDouble> r(10, 0.0); // store det in 10th entry
461 
462  r[9] = 1.0 / (J(0, 0) * (J(1, 1) * J(2, 2) - J(2, 1) * J(1, 2)) -
463  J(0, 1) * (J(1, 0) * J(2, 2) - J(2, 0) * J(1, 2)) +
464  J(0, 2) * (J(1, 0) * J(2, 1) - J(2, 0) * J(1, 1)));
465 
466  r[0] = J(0, 0);
467  r[1] = J(1, 0);
468  r[2] = J(2, 0);
469  r[3] = J(0, 1);
470  r[4] = J(1, 1);
471  r[5] = J(2, 1);
472  r[6] = J(0, 2);
473  r[7] = J(1, 2);
474  r[8] = J(2, 2);
475  mapsStd.push_back(r);
476  }
477 
478  for (int i = 0; i < m_derivUtil->pts; ++i)
479  {
480  NekDouble a1 = 0.5 * (1 - u2[i]);
481  NekDouble a2 = 0.5 * (1 + u2[i]);
482  NekDouble b1 = 0.5 * (1 - v2[i]);
483  NekDouble b2 = 0.5 * (1 + v2[i]);
484  NekDouble c1 = 0.5 * (1 - w2[i]);
485  NekDouble c2 = 0.5 * (1 + w2[i]);
486 
487  DNekMat J(3, 3, 1.0, eFULL);
488 
489  J(0, 0) = -0.5 * b1 * c1 * xyz[0][0] + 0.5 * b1 * c1 * xyz[1][0] +
490  0.5 * b2 * c1 * xyz[2][0] - 0.5 * b2 * c1 * xyz[3][0] -
491  0.5 * b1 * c2 * xyz[5][0] + 0.5 * b1 * c2 * xyz[5][0] +
492  0.5 * b2 * c2 * xyz[6][0] - 0.5 * b2 * c2 * xyz[7][0];
493  J(1, 0) = -0.5 * b1 * c1 * xyz[0][1] + 0.5 * b1 * c1 * xyz[1][1] +
494  0.5 * b2 * c1 * xyz[2][1] - 0.5 * b2 * c1 * xyz[3][1] -
495  0.5 * b1 * c2 * xyz[5][1] + 0.5 * b1 * c2 * xyz[5][1] +
496  0.5 * b2 * c2 * xyz[6][1] - 0.5 * b2 * c2 * xyz[7][1];
497  J(2, 0) = -0.5 * b1 * c1 * xyz[0][2] + 0.5 * b1 * c1 * xyz[1][2] +
498  0.5 * b2 * c1 * xyz[2][2] - 0.5 * b2 * c1 * xyz[3][2] -
499  0.5 * b1 * c2 * xyz[5][2] + 0.5 * b1 * c2 * xyz[5][2] +
500  0.5 * b2 * c2 * xyz[6][2] - 0.5 * b2 * c2 * xyz[7][2];
501 
502  J(0, 1) = -0.5 * a1 * c1 * xyz[0][0] - 0.5 * a2 * c1 * xyz[1][0] +
503  0.5 * a2 * c1 * xyz[2][0] + 0.5 * a1 * c1 * xyz[3][0] -
504  0.5 * a1 * c2 * xyz[5][0] - 0.5 * a2 * c2 * xyz[5][0] +
505  0.5 * a2 * c2 * xyz[6][0] + 0.5 * a1 * c2 * xyz[7][0];
506  J(1, 1) = -0.5 * a1 * c1 * xyz[0][1] - 0.5 * a2 * c1 * xyz[1][1] +
507  0.5 * a2 * c1 * xyz[2][1] + 0.5 * a1 * c1 * xyz[3][1] -
508  0.5 * a1 * c2 * xyz[5][1] - 0.5 * a2 * c2 * xyz[5][1] +
509  0.5 * a2 * c2 * xyz[6][1] + 0.5 * a1 * c2 * xyz[7][1];
510  J(2, 1) = -0.5 * a1 * c1 * xyz[0][2] - 0.5 * a2 * c1 * xyz[1][2] +
511  0.5 * a2 * c1 * xyz[2][2] + 0.5 * a1 * c1 * xyz[3][2] -
512  0.5 * a1 * c2 * xyz[5][2] - 0.5 * a2 * c2 * xyz[5][2] +
513  0.5 * a2 * c2 * xyz[6][2] + 0.5 * a1 * c2 * xyz[7][2];
514 
515  J(0, 0) = -0.5 * b1 * a1 * xyz[0][0] - 0.5 * b1 * a2 * xyz[1][0] -
516  0.5 * b2 * a2 * xyz[2][0] - 0.5 * b2 * a1 * xyz[3][0] +
517  0.5 * b1 * a1 * xyz[5][0] + 0.5 * b1 * a2 * xyz[5][0] +
518  0.5 * b2 * a2 * xyz[6][0] + 0.5 * b2 * a1 * xyz[7][0];
519  J(1, 0) = -0.5 * b1 * a1 * xyz[0][1] - 0.5 * b1 * a2 * xyz[1][1] -
520  0.5 * b2 * a2 * xyz[2][1] - 0.5 * b2 * a1 * xyz[3][1] +
521  0.5 * b1 * a1 * xyz[5][1] + 0.5 * b1 * a2 * xyz[5][1] +
522  0.5 * b2 * a2 * xyz[6][1] + 0.5 * b2 * a1 * xyz[7][1];
523  J(2, 0) = -0.5 * b1 * a1 * xyz[0][2] - 0.5 * b1 * a2 * xyz[1][2] -
524  0.5 * b2 * a2 * xyz[2][2] - 0.5 * b2 * a1 * xyz[3][2] +
525  0.5 * b1 * a1 * xyz[5][2] + 0.5 * b1 * a2 * xyz[5][2] +
526  0.5 * b2 * a2 * xyz[6][2] + 0.5 * b2 * a1 * xyz[7][2];
527 
528  J.Invert();
529 
530  vector<NekDouble> r(10, 0.0); // store det in 10th entry
531 
532  r[9] = 1.0 / (J(0, 0) * (J(1, 1) * J(2, 2) - J(2, 1) * J(1, 2)) -
533  J(0, 1) * (J(1, 0) * J(2, 2) - J(2, 0) * J(1, 2)) +
534  J(0, 2) * (J(1, 0) * J(2, 1) - J(2, 0) * J(1, 1)));
535 
536  r[0] = J(0, 0);
537  r[1] = J(1, 0);
538  r[2] = J(2, 0);
539  r[3] = J(0, 1);
540  r[4] = J(1, 1);
541  r[5] = J(2, 1);
542  r[6] = J(0, 2);
543  r[7] = J(1, 2);
544  r[8] = J(2, 2);
545  maps.push_back(r);
546  }
547  }
548  else
549  {
550  ASSERTL0(false, "not coded");
551  }
552 }
553 
554 void ElUtil::Evaluate()
555 {
556  NekDouble mx2 = -1.0 * numeric_limits<double>::max();
557  NekDouble mn2 = numeric_limits<double>::max();
558 
559  ASSERTL0(nodes.size() == m_derivUtil->ptsStd, "node count wrong");
560  const int nNodes = nodes.size();
561 
562  if (m_dim == 2)
563  {
564  std::vector<NekDouble> X(nNodes), Y(nNodes);
565  for (int j = 0; j < nNodes; j++)
566  {
567  X[j] = *nodes[j][0];
568  Y[j] = *nodes[j][1];
569  }
570 
571  std::vector<NekDouble> x1i(nNodes), y1i(nNodes);
572  std::vector<NekDouble> x2i(nNodes), y2i(nNodes);
573 
574  Blas::Dgemv(
575  'N', nNodes, nNodes, 1.0, m_derivUtil->VdmDStd[0].GetRawPtr(),
576  nNodes, &X[0], 1, 0.0, &x1i[0], 1.0);
577  Blas::Dgemv(
578  'N', nNodes, nNodes, 1.0, m_derivUtil->VdmDStd[0].GetRawPtr(),
579  nNodes, &Y[0], 1, 0.0, &y1i[0], 1.0);
580  Blas::Dgemv(
581  'N', nNodes, nNodes, 1.0, m_derivUtil->VdmDStd[1].GetRawPtr(),
582  nNodes, &X[0], 1, 0.0, &x2i[0], 1.0);
583  Blas::Dgemv(
584  'N', nNodes, nNodes, 1.0, m_derivUtil->VdmDStd[1].GetRawPtr(),
585  nNodes, &Y[0], 1, 0.0, &y2i[0], 1.0);
586 
587  for (int j = 0; j < nNodes; j++)
588  {
589  NekDouble jacDet = x1i[j] * y2i[j] - x2i[j] * y1i[j];
590  mx2 = max(mx2, jacDet / mapsStd[j][9]);
591  mn2 = min(mn2, jacDet / mapsStd[j][9]);
592  }
593  }
594  else if (m_dim == 3)
595  {
596  std::vector<NekDouble> X(nNodes), Y(nNodes), Z(nNodes);
597  for (int j = 0; j < nNodes; j++)
598  {
599  X[j] = *nodes[j][0];
600  Y[j] = *nodes[j][1];
601  Z[j] = *nodes[j][2];
602  }
603  std::vector<NekDouble> x1i2(nNodes), y1i2(nNodes), z1i2(nNodes);
604  std::vector<NekDouble> x2i2(nNodes), y2i2(nNodes), z2i2(nNodes);
605  std::vector<NekDouble> x3i2(nNodes), y3i2(nNodes), z3i2(nNodes);
606 
607  Blas::Dgemv(
608  'N', nNodes, nNodes, 1.0, m_derivUtil->VdmDStd[0].GetRawPtr(),
609  nNodes, &X[0], 1, 0.0, &x1i2[0], 1.0);
610  Blas::Dgemv(
611  'N', nNodes, nNodes, 1.0, m_derivUtil->VdmDStd[0].GetRawPtr(),
612  nNodes, &Y[0], 1, 0.0, &y1i2[0], 1.0);
613  Blas::Dgemv(
614  'N', nNodes, nNodes, 1.0, m_derivUtil->VdmDStd[0].GetRawPtr(),
615  nNodes, &Z[0], 1, 0.0, &z1i2[0], 1.0);
616  Blas::Dgemv(
617  'N', nNodes, nNodes, 1.0, m_derivUtil->VdmDStd[1].GetRawPtr(),
618  nNodes, &X[0], 1, 0.0, &x2i2[0], 1.0);
619  Blas::Dgemv(
620  'N', nNodes, nNodes, 1.0, m_derivUtil->VdmDStd[1].GetRawPtr(),
621  nNodes, &Y[0], 1, 0.0, &y2i2[0], 1.0);
622  Blas::Dgemv(
623  'N', nNodes, nNodes, 1.0, m_derivUtil->VdmDStd[1].GetRawPtr(),
624  nNodes, &Z[0], 1, 0.0, &z2i2[0], 1.0);
625  Blas::Dgemv(
626  'N', nNodes, nNodes, 1.0, m_derivUtil->VdmDStd[2].GetRawPtr(),
627  nNodes, &X[0], 1, 0.0, &x3i2[0], 1.0);
628  Blas::Dgemv(
629  'N', nNodes, nNodes, 1.0, m_derivUtil->VdmDStd[2].GetRawPtr(),
630  nNodes, &Y[0], 1, 0.0, &y3i2[0], 1.0);
631  Blas::Dgemv(
632  'N', nNodes, nNodes, 1.0, m_derivUtil->VdmDStd[2].GetRawPtr(),
633  nNodes, &Z[0], 1, 0.0, &z3i2[0], 1.0);
634 
635  for (int j = 0; j < nNodes; j++)
636  {
637  NekDouble jacDet =
638  x1i2[j] * (y2i2[j] * z3i2[j] - z2i2[j] * y3i2[j]) -
639  x2i2[j] * (y1i2[j] * z3i2[j] - z1i2[j] * y3i2[j]) +
640  x3i2[j] * (y1i2[j] * z2i2[j] - z1i2[j] * y2i2[j]);
641 
642  mx2 = max(mx2, jacDet / mapsStd[j][9]);
643  mn2 = min(mn2, jacDet / mapsStd[j][9]);
644  }
645  }
646 
647  mtx2.lock();
648  if (mn2 < 0)
649  {
650  m_res->startInv++;
651  }
652  m_res->worstJac = min(m_res->worstJac, (mn2 / mx2));
653  mtx2.unlock();
654 
655  m_scaledJac = (mn2 / mx2);
656 }
657 
658 void ElUtil::InitialMinJac()
659 {
660  NekDouble mx = -1.0 * numeric_limits<double>::max();
661  NekDouble mn = numeric_limits<double>::max();
662 
663  ASSERTL0(nodes.size() == m_derivUtil->ptsStd, "node count wrong");
664 
665  if (m_dim == 2)
666  {
667  NekVector<NekDouble> X(nodes.size()), Y(nodes.size());
668  for (int j = 0; j < nodes.size(); j++)
669  {
670  X(j) = *nodes[j][0];
671  Y(j) = *nodes[j][1];
672  }
673 
674  NekVector<NekDouble> x1i2(m_derivUtil->pts), y1i2(m_derivUtil->pts),
675  x2i2(m_derivUtil->pts), y2i2(m_derivUtil->pts);
676 
677  x1i2 = m_derivUtil->VdmD[0] * X;
678  y1i2 = m_derivUtil->VdmD[0] * Y;
679  x2i2 = m_derivUtil->VdmD[1] * X;
680  y2i2 = m_derivUtil->VdmD[1] * Y;
681 
682  for (int j = 0; j < m_derivUtil->pts; j++)
683  {
684  NekDouble jacDet = x1i2(j) * y2i2(j) - x2i2(j) * y1i2(j);
685  jacDet /= maps[j][9];
686  mx = max(mx, jacDet);
687  mn = min(mn, jacDet);
688  }
689  }
690  else if (m_dim == 3)
691  {
692  NekVector<NekDouble> X(nodes.size()), Y(nodes.size()), Z(nodes.size());
693  for (int j = 0; j < nodes.size(); j++)
694  {
695  X(j) = *nodes[j][0];
696  Y(j) = *nodes[j][1];
697  Z(j) = *nodes[j][2];
698  }
699 
700  NekVector<NekDouble> x1i(m_derivUtil->pts), y1i(m_derivUtil->pts),
701  z1i(m_derivUtil->pts), x2i(m_derivUtil->pts), y2i(m_derivUtil->pts),
702  z2i(m_derivUtil->pts), x3i(m_derivUtil->pts), y3i(m_derivUtil->pts),
703  z3i(m_derivUtil->pts);
704 
705  x1i = m_derivUtil->VdmD[0] * X;
706  y1i = m_derivUtil->VdmD[0] * Y;
707  z1i = m_derivUtil->VdmD[0] * Z;
708  x2i = m_derivUtil->VdmD[1] * X;
709  y2i = m_derivUtil->VdmD[1] * Y;
710  z2i = m_derivUtil->VdmD[1] * Z;
711  x3i = m_derivUtil->VdmD[2] * X;
712  y3i = m_derivUtil->VdmD[2] * Y;
713  z3i = m_derivUtil->VdmD[2] * Z;
714 
715  for (int j = 0; j < m_derivUtil->pts; j++)
716  {
717  DNekMat dxdz(3, 3, 1.0, eFULL);
718  dxdz(0, 0) = x1i(j);
719  dxdz(0, 1) = x2i(j);
720  dxdz(0, 2) = x3i(j);
721  dxdz(1, 0) = y1i(j);
722  dxdz(1, 1) = y2i(j);
723  dxdz(1, 2) = y3i(j);
724  dxdz(2, 0) = z1i(j);
725  dxdz(2, 1) = z2i(j);
726  dxdz(2, 2) = z3i(j);
727 
728  NekDouble jacDet =
729  dxdz(0, 0) *
730  (dxdz(1, 1) * dxdz(2, 2) - dxdz(2, 1) * dxdz(1, 2)) -
731  dxdz(0, 1) *
732  (dxdz(1, 0) * dxdz(2, 2) - dxdz(2, 0) * dxdz(1, 2)) +
733  dxdz(0, 2) *
734  (dxdz(1, 0) * dxdz(2, 1) - dxdz(2, 0) * dxdz(1, 1));
735 
736  mx = max(mx, jacDet);
737  mn = min(mn, jacDet);
738  }
739  }
740 
741  m_minJac = mn;
742 }
743 
744 ElUtilJob *ElUtil::GetJob()
745 {
746  return new ElUtilJob(this);
747 }
748 }
749 }
Nektar::Utilities::DerivUtilSharedPtr
std::shared_ptr< DerivUtil > DerivUtilSharedPtr
Definition: ElUtil.h:50
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::FieldUtils::MappingIdealToRef
vector< DNekMat > MappingIdealToRef(SpatialDomains::GeometrySharedPtr geom, StdRegions::StdExpansionSharedPtr chi)
Definition: ProcessQualityMetric.cpp:118
ProcessVarOpti.h
std
STL namespace.
Nektar::Utilities::ResidualSharedPtr
std::shared_ptr< Residual > ResidualSharedPtr
Definition: ElUtil.h:53
ManagerAccess.h
Nektar::NekMeshUtils::ElementSharedPtr
std::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:49
Nektar::LibUtilities::PointsManager
PointsManagerT & PointsManager(void)
Definition: ManagerAccess.cpp:42
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:59
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Blas::Dgemv
static void Dgemv(const char &trans, const int &m, const int &n, const double &alpha, const double *a, const int &lda, const double *x, const int &incx, const double &beta, double *y, const int &incy)
BLAS level 2: Matrix vector multiply y = A x where A[m x n].
Definition: Blas.hpp:168
ElUtil.h
Nektar::Utilities::mtx2
std::mutex mtx2
Definition: ElUtil.cpp:48
Nektar::LibUtilities::eNodalPrismElec
3D electrostatically spaced points on a Prism
Definition: PointsType.h:75
Generated on Wed Jan 20 2021 03:33:33 for Nektar++ by   doxygen 1.8.13