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
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.
void setOrdinalTagData(std::vector< std::vector< std::vector< int > > > &tagToOrdinal, std::vector< std::vector< int > > &ordinalToTag, const int *tags, const int basisCard, const int tagSize, const int posScDim, const int posScOrd, const int posDfOrd)
Fills ordinalToTag_ and tagToOrdinal_ by basis-specific tag data.
bool basisTagsAreSet_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
FieldContainer< Scalar > V_
Generalized Vandermonde matrix V_{ij} = phis_i(x_j)
EOperator
Enumeration of primitive operators available in Intrepid. Primitive operators act on reconstructed fu...
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...
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
EPointType
Enumeration of types of point distributions in Intrepid.
virtual void getDofCoords(ArrayScalar &DofCoords) const
implements the dofcoords interface
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. ...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Line cell.
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.