Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
StdPointExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File $Source: /usr/sci/projects/Nektar/cvs/Nektar++/library/LocalRegions/PointExp.cpp,v $
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: Definition of a Point expansion
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <StdRegions/StdPointExp.h>
37 
38 namespace Nektar
39 {
40  namespace StdRegions
41  {
42 
43  StdPointExp::StdPointExp()
44  {
45  }
46 
47 
48  StdPointExp::StdPointExp(const LibUtilities::BasisKey &Ba):
49  StdExpansion(Ba.GetNumModes(), 1, Ba),
50  StdExpansion0D(Ba.GetNumModes(),Ba)
51  {
52  }
53 
54 
55  /** \brief Copy Constructor */
56 
57  StdPointExp::StdPointExp(const StdPointExp &T):
58  StdExpansion(T),
59  StdExpansion0D(T)
60  {
61  }
62 
63 
64 
65  StdPointExp::~StdPointExp()
66  {
67  }
68 
69 
70  LibUtilities::ShapeType StdPointExp::v_DetShapeType() const
71  {
72  return LibUtilities::ePoint;
73  }
74 
75 
76  void StdPointExp::v_GetCoords(
77  Array<OneD, NekDouble> &coords_0,
78  Array<OneD, NekDouble> &coords_1,
79  Array<OneD, NekDouble> &coords_2)
80  {
81  Blas::Dcopy(GetNumPoints(0),(m_base[0]->GetZ()).get(),
82  1,&coords_0[0],1);
83  }
84 
85 
86  void StdPointExp::v_BwdTrans(
87  const Array<OneD, const NekDouble>& inarray,
88  Array<OneD, NekDouble> &outarray)
89  {
90  int nquad = m_base[0]->GetNumPoints();
91 
92  if(m_base[0]->Collocation())
93  {
94  Vmath::Vcopy(nquad, inarray, 1, outarray, 1);
95  }
96  else
97  {
98 
99 #ifdef NEKTAR_USING_DIRECT_BLAS_CALLS
100 
101  Blas::Dgemv('N',nquad,m_base[0]->GetNumModes(),1.0, (m_base[0]->GetBdata()).get(),
102  nquad,&inarray[0],1,0.0,&outarray[0],1);
103 
104 #else //NEKTAR_USING_DIRECT_BLAS_CALLS
105 
106  NekVector<NekDouble> in(m_ncoeffs,inarray,eWrapper);
107  NekVector<NekDouble> out(nquad,outarray,eWrapper);
108  NekMatrix<NekDouble> B(nquad,m_ncoeffs,m_base[0]->GetBdata(),eWrapper);
109  out = B * in;
110 
111 #endif //NEKTAR_USING_DIRECT_BLAS_CALLS
112 
113  }
114  }
115 
116 
117  void StdPointExp::v_FwdTrans(const Array<OneD, const NekDouble>& inarray,
118  Array<OneD, NekDouble> &outarray)
119  {
120  if(m_base[0]->Collocation())
121  {
122  Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
123  }
124  else
125  {
126  v_IProductWRTBase(inarray,outarray);
127 
128  // get Mass matrix inverse
129  StdMatrixKey masskey(eInvMass,v_DetShapeType(),*this);
130  DNekMatSharedPtr matsys = GetStdMatrix(masskey);
131 
132  NekVector<NekDouble> in(m_ncoeffs,outarray,eCopy);
133  NekVector<NekDouble> out(m_ncoeffs,outarray,eWrapper);
134 
135  out = (*matsys)*in;
136  }
137  }
138 
139  void StdPointExp::v_FwdTrans_BndConstrained(
140  const Array<OneD, const NekDouble>& inarray,
141  Array<OneD, NekDouble> &outarray)
142  {
143  if(m_base[0]->Collocation())
144  {
145  Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
146  }
147  else
148  {
149  int nInteriorDofs = m_ncoeffs-2;
150  int offset;
151 
152  switch(m_base[0]->GetBasisType())
153  {
154  case LibUtilities::eGLL_Lagrange:
155  {
156  offset = 1;
157  }
158  break;
159  case LibUtilities::eModified_A:
160  case LibUtilities::eModified_B:
161  {
162  offset = 2;
163  }
164  break;
165  default:
166  ASSERTL0(false,"This type of FwdTrans is not defined for this shapex type");
167  }
168 
169  fill(outarray.get(), outarray.get()+m_ncoeffs, 0.0 );
170 
171  outarray[GetVertexMap(0)] = inarray[0];
172  outarray[GetVertexMap(1)] = inarray[m_base[0]->GetNumPoints()-1];
173 
174  if(m_ncoeffs>2)
175  {
176  // ideally, we would like to have tmp0 to be replaced by
177  // outarray (currently MassMatrixOp does not allow aliasing)
178  Array<OneD, NekDouble> tmp0(m_ncoeffs);
179  Array<OneD, NekDouble> tmp1(m_ncoeffs);
180 
181  StdMatrixKey masskey(eMass,v_DetShapeType(),*this);
182  MassMatrixOp(outarray,tmp0,masskey);
183  v_IProductWRTBase(inarray,tmp1);
184 
185  Vmath::Vsub(m_ncoeffs, tmp1, 1, tmp0, 1, tmp1, 1);
186 
187  // get Mass matrix inverse (only of interior DOF)
188  DNekMatSharedPtr matsys = (m_stdStaticCondMatrixManager[masskey])->GetBlock(1,1);
189 
190  Blas::Dgemv('N',nInteriorDofs,nInteriorDofs,1.0, &(matsys->GetPtr())[0],
191  nInteriorDofs,tmp1.get()+offset,1,0.0,outarray.get()+offset,1);
192  }
193  }
194 
195  }
196 
197 
198  void StdPointExp::v_BwdTrans_SumFac(const Array<OneD, const NekDouble>& inarray,
199  Array<OneD, NekDouble> &outarray)
200  {
201  v_BwdTrans(inarray, outarray);
202  }
203 
204 //Inner product
205  void StdPointExp::v_IProductWRTBase(const Array<OneD, const NekDouble>& base,
206  const Array<OneD, const NekDouble>& inarray,
207  Array<OneD, NekDouble> &outarray,
208  int coll_check)
209  {
210  int nquad = m_base[0]->GetNumPoints();
211  Array<OneD, NekDouble> tmp(nquad);
212  Array<OneD, const NekDouble> z = m_base[0]->GetZ();
213  Array<OneD, const NekDouble> w = m_base[0]->GetW();
214 
215  Vmath::Vmul(nquad, inarray, 1, w, 1, tmp, 1);
216 
217  if(coll_check&&m_base[0]->Collocation())
218  {
219  Vmath::Vcopy(nquad, tmp, 1, outarray, 1);
220  }
221  else
222  {
223  Blas::Dgemv('T',nquad,m_ncoeffs,1.0,base.get(),nquad,
224  &tmp[0],1,0.0,outarray.get(),1);
225  }
226  }
227 
228  void StdPointExp::v_IProductWRTBase(const Array<OneD, const NekDouble>& inarray,
229  Array<OneD, NekDouble> &outarray)
230  {
231  v_IProductWRTBase(m_base[0]->GetBdata(),inarray,outarray,1);
232  }
233 
234  void StdPointExp::v_IProductWRTDerivBase(
235  const int dir,
236  const Array<OneD, const NekDouble>& inarray,
237  Array<OneD, NekDouble> & outarray)
238  {
239  ASSERTL1(dir >= 0 && dir < 1,"input dir is out of range");
240  v_IProductWRTBase(m_base[0]->GetDbdata(),inarray,outarray,1);
241  }
242 
243  void StdPointExp::v_IProductWRTBase_SumFac(
244  const Array<OneD, const NekDouble>& inarray,
245  Array<OneD, NekDouble> &outarray)
246  {
247  v_IProductWRTBase(m_base[0]->GetBdata(),inarray,outarray,1);
248  }
249 
250 
251  DNekMatSharedPtr StdPointExp::v_GenMatrix(const StdMatrixKey &mkey)
252  {
253  DNekMatSharedPtr Mat;
254  MatrixType mattype;
255 
256  switch(mattype = mkey.GetMatrixType())
257  {
258  case eFwdTrans:
259  {
260  Mat = MemoryManager<DNekMat>::AllocateSharedPtr(m_ncoeffs,m_ncoeffs);
261  StdMatrixKey iprodkey(eIProductWRTBase,v_DetShapeType(),*this);
262  DNekMat &Iprod = *GetStdMatrix(iprodkey);
263  StdMatrixKey imasskey(eInvMass,v_DetShapeType(),*this);
264  DNekMat &Imass = *GetStdMatrix(imasskey);
265 
266  (*Mat) = Imass*Iprod;
267 
268  }
269  break;
270  default:
271  {
272  Mat = StdExpansion::CreateGeneralMatrix(mkey);
273  }
274  break;
275  }
276 
277  return Mat;
278  }
279 
280 
281  DNekMatSharedPtr StdPointExp::v_CreateStdMatrix(const StdMatrixKey &mkey)
282  {
283  return v_GenMatrix(mkey);
284  }
285 
286 
287  }
288 }
Generated on Sun Mar 15 2015 00:11:35 for Nektar++ by   doxygen 1.8.1.2