Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
NodalTetElec.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File NodalTetElec.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: 3D Nodal Tet Electrostatic Point Definitions
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <LibUtilities/BasicUtils/ErrorUtil.hpp>
37 #include <LibUtilities/Foundations/NodalTetElec.h>
38 #include <LibUtilities/Foundations/NodalTetElecData.h>
39 #include <LibUtilities/Foundations/NodalUtil.h>
40 #include <LibUtilities/LinearAlgebra/NekVector.hpp>
41 #include <LibUtilities/LinearAlgebra/NekMatrix.hpp>
42 
43 namespace Nektar
44 {
45  namespace LibUtilities
46  {
47 
48  // ////////////////////////////////////////////////////////
49  // Coordinate the nodal tetrahedron electrostatic points
50 
51  void NodalTetElec::CalculatePoints()
52  {
53  // Allocate the storage for points
54  Points<NekDouble>::CalculatePoints();
55 
56  int index=0,isum=0;
57  const int offset = 5; //offset to match Datafile
58  NekDouble b,c,d;
59  unsigned int numPoints = GetNumPoints();
60 
61  // initialize values
62  for(unsigned int i=0; i < numPoints-2; ++i)
63  {
64  index += NodalTetElecNPTS[i];
65  }
66 
67  for(unsigned int i=0; i < NodalTetElecNPTS[numPoints-2]; ++i, ++index)
68  {
69  // 1 Point Symmetry: aaaa
70  if(int(NodalTetElecData[index][0]))
71  {
72  b = NodalTetElecData[index][6];
73  c = NodalTetElecData[index][7];
74  d = NodalTetElecData[index][8];
75 
76  m_points[0][isum] = 2.0*b - 1.0;
77  m_points[1][isum] = 2.0*c - 1.0;
78  m_points[2][isum] = 2.0*d - 1.0;
79  isum++;
80  continue;
81  }//end symmetry 1
82 
83 
84  // 4 Point symmetry: aaab or abbb
85  if(int(NodalTetElecData[index][1]))
86  {
87  for(unsigned int j=0; j < 4; ++j)
88  {
89  b = NodalTetElecData[index][offset + perm4_3d[j][1]];
90  c = NodalTetElecData[index][offset + perm4_3d[j][2]];
91  d = NodalTetElecData[index][offset + perm4_3d[j][3]];
92 
93  m_points[0][isum] = 2.0*b - 1.0;
94  m_points[1][isum] = 2.0*c - 1.0;
95  m_points[2][isum] = 2.0*d - 1.0;
96  isum++;
97  }//end j
98  continue;
99  }//end symmetry 4
100 
101 
102  // 6 Point symmetry: aabb
103  if(int(NodalTetElecData[index][2]))
104  {
105  for(unsigned int j=0; j < 6; ++j)
106  {
107  b = NodalTetElecData[index][offset + perm6_3d[j][1]];
108  c = NodalTetElecData[index][offset + perm6_3d[j][2]];
109  d = NodalTetElecData[index][offset + perm6_3d[j][3]];
110 
111  m_points[0][isum] = 2.0*b - 1.0;
112  m_points[1][isum] = 2.0*c - 1.0;
113  m_points[2][isum] = 2.0*d - 1.0;
114  isum++;
115  }//end j
116  continue;
117  }//end symmetry6
118 
119 
120  // 12 Point symmetry: case aabc
121  if(int(NodalTetElecData[index][3]) == 1)
122  {
123  for(unsigned int j=0; j < 12; ++j)
124  {
125  b = NodalTetElecData[index][offset + perm12A_3d[j][1]];
126  c = NodalTetElecData[index][offset + perm12A_3d[j][2]];
127  d = NodalTetElecData[index][offset + perm12A_3d[j][3]];
128 
129  m_points[0][isum] = 2.0*b - 1.0;
130  m_points[1][isum] = 2.0*c - 1.0;
131  m_points[2][isum] = 2.0*d - 1.0;
132  isum++;
133  }//end j
134  continue;
135  }//end symmetry 12 aabc
136 
137 
138  // 12 Point symmetry: case abcc
139  if(int(NodalTetElecData[index][3]) == 2)
140  {
141  for(unsigned int j=0; j < 12; ++j)
142  {
143  b = NodalTetElecData[index][offset + perm12B_3d[j][1]];
144  c = NodalTetElecData[index][offset + perm12B_3d[j][2]];
145  d = NodalTetElecData[index][offset + perm12B_3d[j][3]];
146 
147  m_points[0][isum] = 2.0*b - 1.0;
148  m_points[1][isum] = 2.0*c - 1.0;
149  m_points[2][isum] = 2.0*d - 1.0;
150  isum++;
151  }//end j
152  continue;
153  }//end symmetry 12 abcc
154 
155 
156  // 12 Point symmetry: case abbc
157  if(int(NodalTetElecData[index][3]) == 3)
158  {
159  for(unsigned int j=0; j < 12; ++j)
160  {
161  b = NodalTetElecData[index][offset + perm12C_3d[j][1]];
162  c = NodalTetElecData[index][offset + perm12C_3d[j][2]];
163  d = NodalTetElecData[index][offset + perm12C_3d[j][3]];
164 
165  m_points[0][isum] = 2.0*b - 1.0;
166  m_points[1][isum] = 2.0*c - 1.0;
167  m_points[2][isum] = 2.0*d - 1.0;
168  isum++;
169  }//end j
170  continue;
171  }//end symmetry 12 abbc
172 
173  // 24 Point symmetry: case abcd
174  if(int(NodalTetElecData[index][4]))
175  {
176  for(unsigned int j=0; j < 24; ++j)
177  {
178  b = NodalTetElecData[index][offset + perm24_3d[j][1]];
179  c = NodalTetElecData[index][offset + perm24_3d[j][2]];
180  d = NodalTetElecData[index][offset + perm24_3d[j][3]];
181 
182  m_points[0][isum] = 2.0*b - 1.0;
183  m_points[1][isum] = 2.0*c - 1.0;
184  m_points[2][isum] = 2.0*d - 1.0;
185  isum++;
186  }//end j
187  continue;
188  }//end symmetry24abcd
189 
190 
191  }//end npts
192 
193  NodalPointReorder3d();
194 
195  ASSERTL1((static_cast<unsigned int>(isum)==m_pointsKey.GetTotNumPoints()),"sum not equal to npts");
196  }
197 
198  void NodalTetElec::CalculateWeights()
199  {
200  // Allocate the storage for points
201  PointsBaseType::CalculateWeights();
202 
203  typedef DataType T;
204 
205  // Solve the Vandermonde system of integrals for the weight vector
206  NekVector<T> w = MakeTetWeights(NekVector<T>(m_points[0]), NekVector<T>(m_points[1]), NekVector<T>(m_points[2]));
207 
208  m_weights = Array<OneD,T>( w.GetRows(), w.GetPtr() );
209  }
210 
211  // ////////////////////////////////////////
212  // CalculateInterpMatrix()
213  void NodalTetElec::CalculateInterpMatrix(const Array<OneD, const NekDouble>& xia, const Array<OneD, const NekDouble>& yia,
214  const Array<OneD, const NekDouble>& zia, Array<OneD, NekDouble>& interp)
215 
216  {
217  NekVector<NekDouble> x( m_points[0] );
218  NekVector<NekDouble> y( m_points[1] );
219  NekVector<NekDouble> z( m_points[2] );
220  NekVector<NekDouble> xi( xia );
221  NekVector<NekDouble> yi( yia );
222  NekVector<NekDouble> zi( zia );
223  NekMatrix<NekDouble> interMat = GetTetInterpolationMatrix(x, y, z, xi, yi, zi);
224 
225  int rows = xi.GetRows(), cols = GetTotNumPoints();
226  for( int i = 0; i < rows; ++i ) {
227  for( int j = 0; j < cols; ++j ) {
228  interp[j + i*cols] = interMat(i,j);
229  }
230  }
231  }
232 
233  void NodalTetElec::CalculateDerivMatrix()
234  {
235  // Allocate the derivative matrix.
236  PointsBaseType::CalculateDerivMatrix();
237 
238  NekVector<NekDouble> x( m_points[0] );
239  NekVector<NekDouble> y( m_points[1] );
240  NekVector<NekDouble> z( m_points[2] );
241  NekVector<NekDouble> xi = x;
242  NekVector<NekDouble> yi = y;
243  NekVector<NekDouble> zi = z;
244 
245  *m_derivmatrix[0] = *GetTetXDerivativeMatrix(x,y,z,xi,yi,zi);
246 
247  *m_derivmatrix[1] = *GetTetYDerivativeMatrix(x,y,z,xi,yi,zi);
248 
249  *m_derivmatrix[2] = *GetTetZDerivativeMatrix(x,y,z,xi,yi,zi);
250  }
251 
252  boost::shared_ptr<PointsBaseType> NodalTetElec::Create(const PointsKey &key)
253  {
254  boost::shared_ptr<PointsBaseType> returnval(MemoryManager<NodalTetElec>::AllocateSharedPtr(key));
255  returnval->Initialize();
256  return returnval;
257  }
258 
259  void NodalTetElec::NodalPointReorder3d()
260  {
261  }
262 
263  } // end of namespace stdregion
264 } // end of namespace stdregion
265 
266 
Generated on Sun Mar 15 2015 00:11:33 for Nektar++ by   doxygen 1.8.1.2