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
51namespace 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>
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
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.
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
Basis_HGRAD_LINE_Cn_FEM(int order, const ArrayScalar &pts)
Constructor.
FieldContainer< Scalar > latticePts_
Holds the points defining the Lagrange basis.
FieldContainer< Scalar > Vinv_
inverse of Generalized Vandermonde matrix, whose columns store the expansion coefficients of the noda...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Line cell.
FieldContainer< Scalar > V_
Generalized Vandermonde matrix V_{ij} = phis_i(x_j)
virtual void getDofCoords(ArrayScalar &DofCoords) const
implements the dofcoords interface
Basis_HGRAD_LINE_Cn_FEM_JACOBI< Scalar, FieldContainer< Scalar > > Phis_
orthogonal basis
bool basisTagsAreSet_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
int basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom.
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
EBasis basisType_
Type of the basis.
int basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package http://trilinos....