Intrepid
Intrepid_HGRAD_LINE_Cn_FEMDef.hpp
1 #ifndef INTREPID_HGRAD_LINE_CN_FEMDEF_HPP
2 #define INTREPID_HGRAD_LINE_CN_FEMDEF_HPP
3 // @HEADER
4 // ************************************************************************
5 //
6 // Intrepid Package
7 // Copyright (2007) Sandia Corporation
8 //
9 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10 // license for use of this work by or on behalf of the U.S. Government.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Pavel Bochev (pbboche@sandia.gov)
40 // Denis Ridzal (dridzal@sandia.gov), or
41 // Kara Peterson (kjpeter@sandia.gov)
42 //
43 // ************************************************************************
44 // @HEADER
45 
51 namespace Intrepid {
52 
53  template<class Scalar, class ArrayScalar>
55  const ArrayScalar &pts ):
56  latticePts_( n+1 , 1 ),
57  Phis_( n ),
58  V_(n+1,n+1),
59  Vinv_(n+1,n+1)
60  {
61  const int N = n+1;
62  this -> basisCardinality_ = N;
63  this -> basisDegree_ = n;
64  this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >() );
65  this -> basisType_ = BASIS_FEM_FIAT;
66  this -> basisCoordinates_ = COORDINATES_CARTESIAN;
67  this -> basisTagsAreSet_ = false;
68 
69 
70  // check validity of points
71  for (int i=0;i<n;i++) {
72  TEUCHOS_TEST_FOR_EXCEPTION( pts(i,0) >= pts(i+1,0) ,
73  std::runtime_error ,
74  "Intrepid::Basis_HGRAD_LINE_Cn_FEM Illegal points given to constructor" );
75  }
76 
77  // copy points int latticePts, correcting endpoints if needed
78  if (std::abs(pts(0,0)+1.0) < INTREPID_TOL) {
79  latticePts_(0,0) = -1.0;
80  }
81  else {
82  latticePts_(0,0) = pts(0,0);
83  }
84  for (int i=1;i<n;i++) {
85  latticePts_(i,0) = pts(i,0);
86  }
87  if (std::abs(pts(n,0)-1.0) < INTREPID_TOL) {
88  latticePts_(n,0) = 1.0;
89  }
90  else {
91  latticePts_(n,0) = pts(n,0);
92  }
93 
94  // form Vandermonde matrix. Actually, this is the transpose of the VDM,
95  // so we transpose on copy below.
96 
97  Phis_.getValues( V_ , latticePts_ , OPERATOR_VALUE );
98 
99  // now I need to copy V into a Teuchos array to do the inversion
100  Teuchos::SerialDenseMatrix<int,Scalar> Vsdm(N,N);
101  for (int i=0;i<N;i++) {
102  for (int j=0;j<N;j++) {
103  Vsdm(i,j) = V_(i,j);
104  }
105  }
106 
107  // invert the matrix
108  Teuchos::SerialDenseSolver<int,Scalar> solver;
109  solver.setMatrix( rcp( &Vsdm , false ) );
110  solver.invert( );
111 
112  // now I need to copy the inverse into Vinv
113  for (int i=0;i<N;i++) {
114  for (int j=0;j<N;j++) {
115  Vinv_(i,j) = Vsdm(j,i);
116  }
117  }
118 
119  }
120 
121  template<class Scalar, class ArrayScalar>
123  const EPointType &pointType ):
124  latticePts_( n+1 , 1 ),
125  Phis_( n ),
126  V_(n+1,n+1),
127  Vinv_(n+1,n+1)
128  {
129  const int N = n+1;
130  this -> basisCardinality_ = N;
131  this -> basisDegree_ = n;
132  this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >() );
133  this -> basisType_ = BASIS_FEM_FIAT;
134  this -> basisCoordinates_ = COORDINATES_CARTESIAN;
135  this -> basisTagsAreSet_ = false;
136 
137  switch(pointType) {
138  case POINTTYPE_EQUISPACED:
139  PointTools::getLattice<Scalar,FieldContainer<Scalar> >( latticePts_ , this->basisCellTopology_ , n , 0 , POINTTYPE_EQUISPACED );
140  break;
141  case POINTTYPE_SPECTRAL:
142  PointTools::getLattice<Scalar,FieldContainer<Scalar> >( latticePts_ , this->basisCellTopology_ , n , 0 , POINTTYPE_WARPBLEND );
143  break;
144  case POINTTYPE_SPECTRAL_OPEN:
145  PointTools::getGaussPoints<Scalar,FieldContainer<Scalar> >( latticePts_ , n );
146  break;
147  default:
148  TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument , "Basis_HGRAD_LINE_Cn_FEM:: invalid point type" );
149  break;
150  }
151 
152  // form Vandermonde matrix. Actually, this is the transpose of the VDM,
153  // so we transpose on copy below.
154 
155  Phis_.getValues( V_ , latticePts_ , OPERATOR_VALUE );
156 
157  // now I need to copy V into a Teuchos array to do the inversion
158  Teuchos::SerialDenseMatrix<int,Scalar> Vsdm(N,N);
159  for (int i=0;i<N;i++) {
160  for (int j=0;j<N;j++) {
161  Vsdm(i,j) = V_(i,j);
162  }
163  }
164 
165  // invert the matrix
166  Teuchos::SerialDenseSolver<int,Scalar> solver;
167  solver.setMatrix( rcp( &Vsdm , false ) );
168  solver.invert( );
169 
170  // now I need to copy the inverse into Vinv
171  for (int i=0;i<N;i++) {
172  for (int j=0;j<N;j++) {
173  Vinv_(i,j) = Vsdm(j,i);
174  }
175  }
176  }
177 
178 
179  template<class Scalar, class ArrayScalar>
181 
182  // Basis-dependent initializations
183  int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
184  int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
185  int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
186  int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
187 
188  // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
189 
190  int *tags = new int[ tagSize * this->getCardinality() ];
191 
192  int internal_dof;
193  int edge_dof;
194 
195  const int n = this->getDegree();
196 
197  // now we check the points for association
198  if (latticePts_(0,0) == -1.0) {
199  tags[0] = 0;
200  tags[1] = 0;
201  tags[2] = 0;
202  tags[3] = 1;
203  edge_dof = 1;
204  internal_dof = n-1;
205  }
206  else {
207  tags[0] = 1;
208  tags[1] = 0;
209  tags[2] = 0;
210  tags[3] = n+1;
211  edge_dof = 0;
212  internal_dof = n+1;
213  }
214  for (int i=1;i<n;i++) {
215  tags[4*i] = 1;
216  tags[4*i+1] = 0;
217  tags[4*i+2] = -edge_dof + i;
218  tags[4*i+3] = internal_dof;
219  }
220  if (latticePts_(n,0) == 1.0) {
221  tags[4*n] = 0;
222  tags[4*n+1] = 1;
223  tags[4*n+2] = 0;
224  tags[4*n+3] = 1;
225  }
226  else {
227  tags[4*n] = 1;
228  tags[4*n+1] = 0;
229  tags[4*n+2] = n;
230  tags[4*n+3] = n;
231  }
232 
233  Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
234  this -> ordinalToTag_,
235  tags,
236  this -> basisCardinality_,
237  tagSize,
238  posScDim,
239  posScOrd,
240  posDfOrd);
241 
242  delete []tags;
243 
244  }
245 
246 
247 
248  template<class Scalar, class ArrayScalar>
250  const ArrayScalar & inputPoints,
251  const EOperator operatorType) const {
252 
253  // Verify arguments
254 #ifdef HAVE_INTREPID_DEBUG
255  Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
256  inputPoints,
257  operatorType,
258  this -> getBaseCellTopology(),
259  this -> getCardinality() );
260 #endif
261  const int numPts = inputPoints.dimension(0);
262  const int numBf = this->getCardinality();
263 
264  try {
265  switch (operatorType) {
266  case OPERATOR_VALUE:
267  {
268  FieldContainer<Scalar> phisCur( numBf , numPts );
269  Phis_.getValues( phisCur , inputPoints , operatorType );
270  for (int i=0;i<outputValues.dimension(0);i++) {
271  for (int j=0;j<outputValues.dimension(1);j++) {
272  outputValues(i,j) = 0.0;
273  for (int k=0;k<this->getCardinality();k++) {
274  outputValues(i,j) += this->Vinv_(k,i) * phisCur(k,j);
275  }
276  }
277  }
278  }
279  break;
280  case OPERATOR_GRAD:
281  case OPERATOR_D1:
282  case OPERATOR_D2:
283  case OPERATOR_D3:
284  case OPERATOR_D4:
285  case OPERATOR_D5:
286  case OPERATOR_D6:
287  case OPERATOR_D7:
288  case OPERATOR_D8:
289  case OPERATOR_D9:
290  case OPERATOR_D10:
291  {
292  const int dkcard =
293  (operatorType == OPERATOR_GRAD)? getDkCardinality(OPERATOR_D1,1): getDkCardinality(operatorType,1);
294 
295  FieldContainer<Scalar> phisCur( numBf , numPts , dkcard );
296  Phis_.getValues( phisCur , inputPoints , operatorType );
297 
298  for (int i=0;i<outputValues.dimension(0);i++) {
299  for (int j=0;j<outputValues.dimension(1);j++) {
300  for (int k=0;k<outputValues.dimension(2);k++) {
301  outputValues(i,j,k) = 0.0;
302  for (int l=0;l<this->getCardinality();l++) {
303  outputValues(i,j,k) += this->Vinv_(l,i) * phisCur(l,j,k);
304  }
305  }
306  }
307  }
308  }
309  break;
310  default:
311  TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
312  ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM): Operator type not implemented" );
313  break;
314  }
315  }
316  catch (std::invalid_argument &exception){
317  TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
318  ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM): Operator failed");
319  }
320 
321  }
322 
323 
324 
325  template<class Scalar, class ArrayScalar>
327  const ArrayScalar & inputPoints,
328  const ArrayScalar & cellVertices,
329  const EOperator operatorType) const {
330  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
331  ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM): FEM Basis calling an FVD member function");
332  }
333 
334 
335  template<class Scalar, class ArrayScalar>
336  void Basis_HGRAD_LINE_Cn_FEM<Scalar,ArrayScalar>::getDofCoords( ArrayScalar & dofCoords ) const
337  {
338  for (int i=0;i<latticePts_.dimension(0);i++)
339  {
340  for (int j=0;j<latticePts_.dimension(1);j++)
341  {
342  dofCoords(i,j) = latticePts_(i,j);
343  }
344  }
345  return;
346  }
347 
348 }// namespace Intrepid
349 #endif
350 
351 #if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
352 #ifdef __GNUC__
353 #warning "The Intrepid package is deprecated"
354 #endif
355 #endif
356 
FieldContainer< Scalar > latticePts_
Holds the points defining the Lagrange basis.
Basis_HGRAD_LINE_Cn_FEM(int order, const ArrayScalar &pts)
Constructor.
EBasis basisType_
Type of the basis.
bool basisTagsAreSet_
&quot;true&quot; if tagToOrdinal_ and ordinalToTag_ have been initialized
FieldContainer< Scalar > V_
Generalized Vandermonde matrix V_{ij} = phis_i(x_j)
Basis_HGRAD_LINE_Cn_FEM_JACOBI< Scalar, FieldContainer< Scalar > > Phis_
orthogonal basis
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package http://trilinos...
virtual void getDofCoords(ArrayScalar &DofCoords) const
implements the dofcoords interface
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Line cell.
int basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
int basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom. ...
FieldContainer< Scalar > Vinv_
inverse of Generalized Vandermonde matrix, whose columns store the expansion coefficients of the noda...
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.