Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Calculate jacobians of elements.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include "ElUtil.h"
37 #include "ProcessVarOpti.h"
38 
39 #include <LibUtilities/Foundations/ManagerAccess.h>
40 
41 using namespace std;
42 
43 namespace Nektar
44 {
45 namespace Utilities
46 {
47 
48 boost::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<Array<OneD, NekDouble> > xyz(4);
97  vector<NodeSharedPtr> ns = m_el->GetVertexList();
98  for (int i = 0; i < 4; i++)
99  {
100  Array<OneD, 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  Array<OneD, 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  Array<OneD, 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  Array<OneD, 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  Array<OneD, 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<Array<OneD, NekDouble> > xyz(6);
268  vector<NodeSharedPtr> ns = m_el->GetVertexList();
269  for (int i = 0; i < 6; i++)
270  {
271  Array<OneD, 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  Array<OneD, 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  Array<OneD, 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<Array<OneD, NekDouble> > xyz(8);
398  vector<NodeSharedPtr> ns = m_el->GetVertexList();
399  for (int i = 0; i < 8; i++)
400  {
401  Array<OneD, 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  Array<OneD, 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  Array<OneD, 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 
561  if (m_dim == 2)
562  {
563  NekVector<NekDouble> X(nodes.size()), Y(nodes.size());
564  for (int j = 0; j < nodes.size(); j++)
565  {
566  X(j) = *nodes[j][0];
567  Y(j) = *nodes[j][1];
568  }
569 
570  NekVector<NekDouble> x1i(nodes.size()), y1i(nodes.size()),
571  x2i(nodes.size()), y2i(nodes.size());
572 
573  x1i = m_derivUtil->VdmDStd[0] * X;
574  y1i = m_derivUtil->VdmDStd[0] * Y;
575  x2i = m_derivUtil->VdmDStd[1] * X;
576  y2i = m_derivUtil->VdmDStd[1] * Y;
577 
578  for (int j = 0; j < nodes.size(); j++)
579  {
580  NekDouble jacDet = x1i(j) * y2i(j) - x2i(j) * y1i(j);
581  mx2 = max(mx2, jacDet / mapsStd[j][9]);
582  mn2 = min(mn2, jacDet / mapsStd[j][9]);
583  }
584  }
585  else if (m_dim == 3)
586  {
587  NekVector<NekDouble> X(nodes.size()), Y(nodes.size()), Z(nodes.size());
588  for (int j = 0; j < nodes.size(); j++)
589  {
590  X(j) = *nodes[j][0];
591  Y(j) = *nodes[j][1];
592  Z(j) = *nodes[j][2];
593  }
594  NekVector<NekDouble> x1i2(nodes.size()), y1i2(nodes.size()),
595  z1i2(nodes.size()), x2i2(nodes.size()), y2i2(nodes.size()),
596  z2i2(nodes.size()), x3i2(nodes.size()), y3i2(nodes.size()),
597  z3i2(nodes.size());
598 
599  x1i2 = m_derivUtil->VdmDStd[0] * X;
600  y1i2 = m_derivUtil->VdmDStd[0] * Y;
601  z1i2 = m_derivUtil->VdmDStd[0] * Z;
602  x2i2 = m_derivUtil->VdmDStd[1] * X;
603  y2i2 = m_derivUtil->VdmDStd[1] * Y;
604  z2i2 = m_derivUtil->VdmDStd[1] * Z;
605  x3i2 = m_derivUtil->VdmDStd[2] * X;
606  y3i2 = m_derivUtil->VdmDStd[2] * Y;
607  z3i2 = m_derivUtil->VdmDStd[2] * Z;
608 
609  for (int j = 0; j < nodes.size(); j++)
610  {
611  DNekMat dxdz(3, 3, 1.0, eFULL);
612  dxdz(0, 0) = x1i2(j);
613  dxdz(0, 1) = x2i2(j);
614  dxdz(0, 2) = x3i2(j);
615  dxdz(1, 0) = y1i2(j);
616  dxdz(1, 1) = y2i2(j);
617  dxdz(1, 2) = y3i2(j);
618  dxdz(2, 0) = z1i2(j);
619  dxdz(2, 1) = z2i2(j);
620  dxdz(2, 2) = z3i2(j);
621 
622  NekDouble jacDet =
623  dxdz(0, 0) *
624  (dxdz(1, 1) * dxdz(2, 2) - dxdz(2, 1) * dxdz(1, 2)) -
625  dxdz(0, 1) *
626  (dxdz(1, 0) * dxdz(2, 2) - dxdz(2, 0) * dxdz(1, 2)) +
627  dxdz(0, 2) *
628  (dxdz(1, 0) * dxdz(2, 1) - dxdz(2, 0) * dxdz(1, 1));
629 
630  mx2 = max(mx2, jacDet / mapsStd[j][9]);
631  mn2 = min(mn2, jacDet / mapsStd[j][9]);
632  }
633  }
634 
635  mtx2.lock();
636  if (mn2 < 0)
637  {
638  m_res->startInv++;
639  }
640  m_res->worstJac = min(m_res->worstJac, (mn2 / mx2));
641  mtx2.unlock();
642 
643  m_scaledJac = (mn2 / mx2);
644 }
645 
646 void ElUtil::InitialMinJac()
647 {
648  NekDouble mx = -1.0 * numeric_limits<double>::max();
649  NekDouble mn = numeric_limits<double>::max();
650 
651  ASSERTL0(nodes.size() == m_derivUtil->ptsStd, "node count wrong");
652 
653  if (m_dim == 2)
654  {
655  NekVector<NekDouble> X(nodes.size()), Y(nodes.size());
656  for (int j = 0; j < nodes.size(); j++)
657  {
658  X(j) = *nodes[j][0];
659  Y(j) = *nodes[j][1];
660  }
661 
662  NekVector<NekDouble> x1i2(m_derivUtil->pts), y1i2(m_derivUtil->pts),
663  x2i2(m_derivUtil->pts), y2i2(m_derivUtil->pts);
664 
665  x1i2 = m_derivUtil->VdmD[0] * X;
666  y1i2 = m_derivUtil->VdmD[0] * Y;
667  x2i2 = m_derivUtil->VdmD[1] * X;
668  y2i2 = m_derivUtil->VdmD[1] * Y;
669 
670  for (int j = 0; j < m_derivUtil->pts; j++)
671  {
672  NekDouble jacDet = x1i2(j) * y2i2(j) - x2i2(j) * y1i2(j);
673  jacDet /= maps[j][9];
674  mx = max(mx, jacDet);
675  mn = min(mn, jacDet);
676  }
677  }
678  else if (m_dim == 3)
679  {
680  NekVector<NekDouble> X(nodes.size()), Y(nodes.size()), Z(nodes.size());
681  for (int j = 0; j < nodes.size(); j++)
682  {
683  X(j) = *nodes[j][0];
684  Y(j) = *nodes[j][1];
685  Z(j) = *nodes[j][2];
686  }
687 
688  NekVector<NekDouble> x1i(m_derivUtil->pts), y1i(m_derivUtil->pts),
689  z1i(m_derivUtil->pts), x2i(m_derivUtil->pts), y2i(m_derivUtil->pts),
690  z2i(m_derivUtil->pts), x3i(m_derivUtil->pts), y3i(m_derivUtil->pts),
691  z3i(m_derivUtil->pts);
692 
693  x1i = m_derivUtil->VdmD[0] * X;
694  y1i = m_derivUtil->VdmD[0] * Y;
695  z1i = m_derivUtil->VdmD[0] * Z;
696  x2i = m_derivUtil->VdmD[1] * X;
697  y2i = m_derivUtil->VdmD[1] * Y;
698  z2i = m_derivUtil->VdmD[1] * Z;
699  x3i = m_derivUtil->VdmD[2] * X;
700  y3i = m_derivUtil->VdmD[2] * Y;
701  z3i = m_derivUtil->VdmD[2] * Z;
702 
703  for (int j = 0; j < m_derivUtil->pts; j++)
704  {
705  DNekMat dxdz(3, 3, 1.0, eFULL);
706  dxdz(0, 0) = x1i(j);
707  dxdz(0, 1) = x2i(j);
708  dxdz(0, 2) = x3i(j);
709  dxdz(1, 0) = y1i(j);
710  dxdz(1, 1) = y2i(j);
711  dxdz(1, 2) = y3i(j);
712  dxdz(2, 0) = z1i(j);
713  dxdz(2, 1) = z2i(j);
714  dxdz(2, 2) = z3i(j);
715 
716  NekDouble jacDet =
717  dxdz(0, 0) *
718  (dxdz(1, 1) * dxdz(2, 2) - dxdz(2, 1) * dxdz(1, 2)) -
719  dxdz(0, 1) *
720  (dxdz(1, 0) * dxdz(2, 2) - dxdz(2, 0) * dxdz(1, 2)) +
721  dxdz(0, 2) *
722  (dxdz(1, 0) * dxdz(2, 1) - dxdz(2, 0) * dxdz(1, 1));
723 
724  mx = max(mx, jacDet);
725  mn = min(mn, jacDet);
726  }
727  }
728 
729  m_minJac = mn;
730 }
731 
732 ElUtilJob *ElUtil::GetJob()
733 {
734  return new ElUtilJob(this);
735 }
736 }
737 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Nektar::FieldUtils::MappingIdealToRef
vector< DNekMat > MappingIdealToRef(SpatialDomains::GeometrySharedPtr geom, StdRegions::StdExpansionSharedPtr chi)
Definition: ProcessQualityMetric.cpp:115
ProcessVarOpti.h
Nektar::Utilities::mtx2
boost::mutex mtx2
Definition: ElUtil.cpp:48
std
STL namespace.
Nektar::Utilities::ResidualSharedPtr
boost::shared_ptr< Residual > ResidualSharedPtr
Definition: ElUtil.h:54
ManagerAccess.h
Nektar::LibUtilities::PointsManager
PointsManagerT & PointsManager(void)
Definition: ManagerAccess.cpp:110
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:58
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
ElUtil.h
Nektar::Utilities::DerivUtilSharedPtr
boost::shared_ptr< DerivUtil > DerivUtilSharedPtr
Definition: ElUtil.h:51
Nektar::NekMeshUtils::ElementSharedPtr
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:49
Vmath::mutex
static boost::mutex mutex
Definition: Vmath.cpp:134
Nektar::LibUtilities::eNodalPrismElec
3D electrostatically spaced points on a Prism
Definition: PointsType.h:76
Generated on Sat Nov 25 2017 16:18:11 for Nektar++ by   doxygen 1.8.8