Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Public Member Functions | Public Attributes | Private Member Functions | Private Attributes | List of all members
Nektar::Utilities::ElUtil Class Reference

#include <ElUtil.h>

Inheritance diagram for Nektar::Utilities::ElUtil:
Inheritance graph
[legend]
Collaboration diagram for Nektar::Utilities::ElUtil:
Collaboration graph
[legend]

Public Member Functions

 ElUtil (ElementSharedPtr e, DerivUtilSharedPtr d, ResidualSharedPtr, int n, int o)
 
ElUtilJobGetJob ()
 
int GetId ()
 
void Evaluate ()
 
void InitialMinJac ()
 
ElementSharedPtr GetEl ()
 
int NodeId (int in)
 
NekDouble GetScaledJac ()
 
NekDoubleGetMinJac ()
 

Public Attributes

std::vector< std::vector
< NekDouble * > > 
nodes
 
std::vector< Array< OneD,
NekDouble > > 
maps
 
std::vector< Array< OneD,
NekDouble > > 
mapsStd
 

Private Member Functions

void MappingIdealToRef ()
 

Private Attributes

ElementSharedPtr m_el
 
int m_dim
 
int m_mode
 
int m_order
 
std::map< int, int > m_idmap
 
NekDouble m_scaledJac
 
NekDouble m_minJac
 
DerivUtilSharedPtr m_derivUtil
 
ResidualSharedPtr m_res
 

Detailed Description

Definition at line 58 of file ElUtil.h.

Constructor & Destructor Documentation

Nektar::Utilities::ElUtil::ElUtil ( ElementSharedPtr  e,
DerivUtilSharedPtr  d,
ResidualSharedPtr  r,
int  n,
int  o 
)

Definition at line 50 of file ElUtil.cpp.

References Nektar::FieldUtils::MappingIdealToRef().

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  }
80 }
std::map< int, int > m_idmap
Definition: ElUtil.h:106
ElementSharedPtr m_el
Definition: ElUtil.h:102
std::vector< std::vector< NekDouble * > > nodes
Definition: ElUtil.h:72
ResidualSharedPtr m_res
Definition: ElUtil.h:112
DerivUtilSharedPtr m_derivUtil
Definition: ElUtil.h:111

Member Function Documentation

void Nektar::Utilities::ElUtil::Evaluate ( )

Definition at line 554 of file ElUtil.cpp.

References ASSERTL0, and Nektar::eFULL.

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 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
boost::mutex mtx2
Definition: ElUtil.cpp:48
std::vector< std::vector< NekDouble * > > nodes
Definition: ElUtil.h:72
ResidualSharedPtr m_res
Definition: ElUtil.h:112
NekMatrix< NekDouble, StandardMatrixTag > DNekMat
Definition: NekTypeDefs.hpp:52
double NekDouble
std::vector< Array< OneD, NekDouble > > mapsStd
Definition: ElUtil.h:73
DerivUtilSharedPtr m_derivUtil
Definition: ElUtil.h:111
ElementSharedPtr Nektar::Utilities::ElUtil::GetEl ( )
inline

Definition at line 78 of file ElUtil.h.

79  {
80  return m_el;
81  }
ElementSharedPtr m_el
Definition: ElUtil.h:102
int Nektar::Utilities::ElUtil::GetId ( )
inline

Definition at line 66 of file ElUtil.h.

67  {
68  return m_el->GetId();
69  }
ElementSharedPtr m_el
Definition: ElUtil.h:102
ElUtilJob * Nektar::Utilities::ElUtil::GetJob ( )

Definition at line 732 of file ElUtil.cpp.

733 {
734  return new ElUtilJob(this);
735 }
NekDouble& Nektar::Utilities::ElUtil::GetMinJac ( )
inline

Definition at line 93 of file ElUtil.h.

94  {
95  return m_minJac;
96  }
NekDouble Nektar::Utilities::ElUtil::GetScaledJac ( )
inline

Definition at line 88 of file ElUtil.h.

89  {
90  return m_scaledJac;
91  }
void Nektar::Utilities::ElUtil::InitialMinJac ( )

Definition at line 646 of file ElUtil.cpp.

References ASSERTL0, and Nektar::eFULL.

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 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
std::vector< std::vector< NekDouble * > > nodes
Definition: ElUtil.h:72
NekMatrix< NekDouble, StandardMatrixTag > DNekMat
Definition: NekTypeDefs.hpp:52
double NekDouble
std::vector< Array< OneD, NekDouble > > maps
Definition: ElUtil.h:73
DerivUtilSharedPtr m_derivUtil
Definition: ElUtil.h:111
void Nektar::Utilities::ElUtil::MappingIdealToRef ( )
private

Definition at line 82 of file ElUtil.cpp.

References ASSERTL0, Nektar::eFULL, Nektar::LibUtilities::eHexahedron, Nektar::LibUtilities::eNodalHexElec, Nektar::LibUtilities::eNodalPrismElec, Nektar::LibUtilities::eNodalPrismSPI, Nektar::LibUtilities::eNodalQuadElec, Nektar::LibUtilities::ePrism, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eTetrahedron, Nektar::LibUtilities::eTriangle, and Nektar::LibUtilities::PointsManager().

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,
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,
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,
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 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
ElementSharedPtr m_el
Definition: ElUtil.h:102
std::vector< std::vector< NekDouble * > > nodes
Definition: ElUtil.h:72
PointsManagerT & PointsManager(void)
NekMatrix< NekDouble, StandardMatrixTag > DNekMat
Definition: NekTypeDefs.hpp:52
double NekDouble
std::vector< Array< OneD, NekDouble > > mapsStd
Definition: ElUtil.h:73
std::vector< Array< OneD, NekDouble > > maps
Definition: ElUtil.h:73
DerivUtilSharedPtr m_derivUtil
Definition: ElUtil.h:111
3D electrostatically spaced points on a Prism
Definition: PointsType.h:76
int Nektar::Utilities::ElUtil::NodeId ( int  in)
inline

Definition at line 83 of file ElUtil.h.

84  {
85  return m_idmap[in];
86  }
std::map< int, int > m_idmap
Definition: ElUtil.h:106

Member Data Documentation

DerivUtilSharedPtr Nektar::Utilities::ElUtil::m_derivUtil
private

Definition at line 111 of file ElUtil.h.

int Nektar::Utilities::ElUtil::m_dim
private

Definition at line 103 of file ElUtil.h.

ElementSharedPtr Nektar::Utilities::ElUtil::m_el
private

Definition at line 102 of file ElUtil.h.

std::map<int,int> Nektar::Utilities::ElUtil::m_idmap
private

Definition at line 106 of file ElUtil.h.

NekDouble Nektar::Utilities::ElUtil::m_minJac
private

Definition at line 109 of file ElUtil.h.

int Nektar::Utilities::ElUtil::m_mode
private

Definition at line 104 of file ElUtil.h.

int Nektar::Utilities::ElUtil::m_order
private

Definition at line 105 of file ElUtil.h.

ResidualSharedPtr Nektar::Utilities::ElUtil::m_res
private

Definition at line 112 of file ElUtil.h.

NekDouble Nektar::Utilities::ElUtil::m_scaledJac
private

Definition at line 108 of file ElUtil.h.

std::vector<Array<OneD, NekDouble> > Nektar::Utilities::ElUtil::maps

Definition at line 73 of file ElUtil.h.

std::vector<Array<OneD, NekDouble> > Nektar::Utilities::ElUtil::mapsStd

Definition at line 73 of file ElUtil.h.

std::vector<std::vector<NekDouble *> > Nektar::Utilities::ElUtil::nodes

Definition at line 72 of file ElUtil.h.