Intrepid
Intrepid_HGRAD_QUAD_Cn_FEMDef.hpp
Go to the documentation of this file.
1#ifndef INTREPID_HGRAD_QUAD_CN_FEMDEF_HPP
2#define INTREPID_HGRAD_QUAD_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
51using Teuchos::Array;
52using Teuchos::RCP;
53
54namespace Intrepid {
55 template<class Scalar, class ArrayScalar>
57 const ArrayScalar &pts_x ,
58 const ArrayScalar &pts_y ):
59 ptsx_( pts_x.dimension(0) , 1 ) ,
60 ptsy_( pts_y.dimension(0) , 1 )
61 {
62 Array<Array<RCP<Basis<Scalar,ArrayScalar> > > > bases(1);
63 bases[0].resize(2);
64 bases[0][0] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM<Scalar,ArrayScalar>( orderx , pts_x ) );
65 bases[0][1] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM<Scalar,ArrayScalar>( ordery , pts_y ) );
66 this->setBases( bases );
67
68 this->basisCardinality_ = (orderx+1)*(ordery+1);
69 if (orderx > ordery) {
70 this->basisDegree_ = orderx;
71 }
72 else {
73 this->basisDegree_ = ordery;
74 }
75 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
76 this -> basisType_ = BASIS_FEM_FIAT;
77 this -> basisCoordinates_ = COORDINATES_CARTESIAN;
78 this -> basisTagsAreSet_ = false;
79
80 for (int i=0;i<pts_x.dimension(0);i++)
81 {
82 ptsx_(i,0) = pts_x(i,0);
83 }
84
85 for (int i=0;i<pts_y.dimension(0);i++)
86 {
87 ptsy_(i,0) = pts_y(i,0);
88 }
89
90 }
91
92 template<class Scalar, class ArrayScalar>
94 const EPointType &pointType ):
95 ptsx_( order+1 , 1 ) ,
96 ptsy_( order+1 , 1 )
97 {
98 Array<Array<RCP<Basis<Scalar,ArrayScalar> > > > bases(1);
99 bases[0].resize(2);
100 bases[0][0] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM<Scalar,ArrayScalar>( order , pointType ) );
101 bases[0][1] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM<Scalar,ArrayScalar>( order , pointType ) );
102 this->setBases( bases );
103
104 this->basisCardinality_ = (order+1)*(order+1);
105 this->basisDegree_ = order;
106 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
107 this -> basisType_ = BASIS_FEM_FIAT;
108 this -> basisCoordinates_ = COORDINATES_CARTESIAN;
109 this -> basisTagsAreSet_ = false;
110
111 // fill up the pt arrays with calls to the lattice
112 EPointType pt = (pointType==POINTTYPE_EQUISPACED)?pointType:POINTTYPE_WARPBLEND;
113 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( ptsx_ ,
114 shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) ,
115 order ,
116 0 ,
117 pt );
118
119 for (int i=0;i<order+1;i++)
120 {
121 ptsy_(i,0) = ptsx_(i,0);
122 }
123 }
124
125 template<class Scalar, class ArrayScalar>
127 {
128 // Basis-dependent initializations
129 int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
130 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
131 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
132 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
133
134 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
135
136 std::vector<int> tags( tagSize * this->getCardinality() );
137
138 // temporarily just put everything on the cell itself
139 for (int i=0;i<this->getCardinality();i++) {
140 tags[tagSize*i] = 2;
141 tags[tagSize*i+1] = 0;
142 tags[tagSize*i+2] = i;
143 tags[tagSize*i+3] = this->getCardinality();
144 }
145
146 Basis<Scalar,ArrayScalar> &xBasis_ = *this->getBases()[0][0];
147 Basis<Scalar,ArrayScalar> &yBasis_ = *this->getBases()[0][1];
148
149 // now let's try to do it "right"
150 // let's get the x and y bases and their dof
151 const std::vector<std::vector<int> >& xdoftags = xBasis_.getAllDofTags();
152 const std::vector<std::vector<int> >& ydoftags = yBasis_.getAllDofTags();
153
154 std::map<int,std::map<int,int> > total_dof_per_entity;
155 std::map<int,std::map<int,int> > current_dof_per_entity;
156
157 for (int i=0;i<4;i++) {
158 total_dof_per_entity[0][i] = 0;
159 current_dof_per_entity[0][i] = 0;
160 }
161 for (int i=0;i<4;i++) {
162 total_dof_per_entity[1][i] = 0;
163 current_dof_per_entity[1][i] = 0;
164 }
165 total_dof_per_entity[2][0] = 0;
166 current_dof_per_entity[2][0] = 0;
167
168 // let's tally the total degrees of freedom on each entity
169 for (int j=0;j<yBasis_.getCardinality();j++) {
170 const int ydim = ydoftags[j][0];
171 const int yent = ydoftags[j][1];
172 for (int i=0;i<xBasis_.getCardinality();i++) {
173 const int xdim = xdoftags[i][0];
174 const int xent = xdoftags[i][1];
175 int dofdim;
176 int dofent;
177 ProductTopology::lineProduct2d( xdim , xent , ydim , yent , dofdim , dofent );
178 total_dof_per_entity[dofdim][dofent] += 1;
179 }
180 }
181
182 int tagcur = 0;
183 for (int j=0;j<yBasis_.getCardinality();j++) {
184 const int ydim = ydoftags[j][0];
185 const int yent = ydoftags[j][1];
186 for (int i=0;i<xBasis_.getCardinality();i++) {
187 const int xdim = xdoftags[i][0];
188 const int xent = xdoftags[i][1];
189 int dofdim;
190 int dofent;
191 ProductTopology::lineProduct2d( xdim , xent , ydim , yent , dofdim , dofent );
192 tags[4*tagcur] = dofdim;
193 tags[4*tagcur+1] = dofent;
194 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
195 current_dof_per_entity[dofdim][dofent]++;
196 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
197 tagcur++;
198 }
199 }
200
201 setOrdinalTagData(this -> tagToOrdinal_,
202 this -> ordinalToTag_,
203 &(tags[0]),
204 this -> basisCardinality_,
205 tagSize,
206 posScDim,
207 posScOrd,
208 posDfOrd);
209
210 }
211
212 template<class Scalar, class ArrayScalar>
214 const ArrayScalar &inputPoints ,
215 const EOperator operatorType ) const
216 {
217#ifdef HAVE_INTREPID_DEBUG
218 getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
219 inputPoints,
220 operatorType,
221 this -> getBaseCellTopology(),
222 this -> getCardinality() );
223#endif
224
225 FieldContainer<Scalar> xInputPoints(inputPoints.dimension(0),1);
226 FieldContainer<Scalar> yInputPoints(inputPoints.dimension(0),1);
227
228 const Basis<Scalar,ArrayScalar> &xBasis_ = *this->bases_[0][0];
229 const Basis<Scalar,ArrayScalar> &yBasis_ = *this->bases_[0][1];
230
231 for (int i=0;i<inputPoints.dimension(0);i++) {
232 xInputPoints(i,0) = inputPoints(i,0);
233 yInputPoints(i,0) = inputPoints(i,1);
234 }
235
236 switch (operatorType) {
237 case OPERATOR_VALUE:
238 {
239 FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0));
240 FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0));
241
242 xBasis_.getValues(xBasisValues,xInputPoints,OPERATOR_VALUE);
243 yBasis_.getValues(yBasisValues,yInputPoints,OPERATOR_VALUE);
244
245 int bfcur = 0;
246 for (int j=0;j<yBasis_.getCardinality();j++) {
247 for (int i=0;i<xBasis_.getCardinality();i++) {
248 for (int k=0;k<inputPoints.dimension(0);k++) {
249 outputValues(bfcur,k) = xBasisValues(i,k) * yBasisValues(j,k);
250 }
251 bfcur++;
252 }
253 }
254 }
255 break;
256 case OPERATOR_GRAD:
257 case OPERATOR_D1:
258 {
259 FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0));
260 FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0));
261 FieldContainer<Scalar> xBasisDerivs(xBasis_.getCardinality(),xInputPoints.dimension(0),1);
262 FieldContainer<Scalar> yBasisDerivs(yBasis_.getCardinality(),yInputPoints.dimension(0),1);
263
264 xBasis_.getValues(xBasisValues,xInputPoints,OPERATOR_VALUE);
265 yBasis_.getValues(yBasisValues,yInputPoints,OPERATOR_VALUE);
266 xBasis_.getValues(xBasisDerivs,xInputPoints,OPERATOR_D1);
267 yBasis_.getValues(yBasisDerivs,yInputPoints,OPERATOR_D1);
268
269 // there are two multiindices: I need the (1,0) and (0,1) derivatives
270 int bfcur = 0;
271
272 for (int j=0;j<yBasis_.getCardinality();j++) {
273 for (int i=0;i<xBasis_.getCardinality();i++) {
274 for (int k=0;k<inputPoints.dimension(0);k++) {
275 outputValues(bfcur,k,0) = xBasisDerivs(i,k,0) * yBasisValues(j,k);
276 outputValues(bfcur,k,1) = xBasisValues(i,k) * yBasisDerivs(j,k,0);
277 }
278 bfcur++;
279 }
280 }
281 }
282 break;
283 case OPERATOR_D2:
284 case OPERATOR_D3:
285 case OPERATOR_D4:
286 case OPERATOR_D5:
287 case OPERATOR_D6:
288 case OPERATOR_D7:
289 case OPERATOR_D8:
290 case OPERATOR_D9:
291 case OPERATOR_D10:
292 {
293 FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0));
294 FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0));
295
296 Teuchos::Array<int> partialMult;
297
298 for (int d=0;d<getDkCardinality(operatorType,2);d++) {
299 getDkMultiplicities( partialMult , d , operatorType , 2 );
300 if (partialMult[0] == 0) {
301 xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0));
302 xBasis_.getValues( xBasisValues , xInputPoints, OPERATOR_VALUE );
303 }
304 else {
305 xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0),1);
306 EOperator xop = (EOperator) ( (int) OPERATOR_D1 + partialMult[0] - 1 );
307 xBasis_.getValues( xBasisValues , xInputPoints, xop );
308 xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0));
309 }
310 if (partialMult[1] == 0) {
311 yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0));
312 yBasis_.getValues( yBasisValues , yInputPoints, OPERATOR_VALUE );
313 }
314 else {
315 yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0),1);
316 EOperator yop = (EOperator) ( (int) OPERATOR_D1 + partialMult[1] - 1 );
317 yBasis_.getValues( yBasisValues , yInputPoints, yop );
318 yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0));
319 }
320
321
322 int bfcur = 0;
323 for (int j=0;j<yBasis_.getCardinality();j++) {
324 for (int i=0;i<xBasis_.getCardinality();i++) {
325 for (int k=0;k<inputPoints.dimension(0);k++) {
326 outputValues(bfcur,k,d) = xBasisValues(i,k) * yBasisValues(j,k);
327 }
328 bfcur++;
329 }
330 }
331 }
332 }
333 break;
334 case OPERATOR_CURL:
335 {
336 FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0));
337 FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0));
338 FieldContainer<Scalar> xBasisDerivs(xBasis_.getCardinality(),xInputPoints.dimension(0),1);
339 FieldContainer<Scalar> yBasisDerivs(yBasis_.getCardinality(),yInputPoints.dimension(0),1);
340
341 xBasis_.getValues(xBasisValues,xInputPoints,OPERATOR_VALUE);
342 yBasis_.getValues(yBasisValues,yInputPoints,OPERATOR_VALUE);
343 xBasis_.getValues(xBasisDerivs,xInputPoints,OPERATOR_D1);
344 yBasis_.getValues(yBasisDerivs,yInputPoints,OPERATOR_D1);
345
346 // there are two multiindices: I need the (1,0) and (0,1) derivatives
347 int bfcur = 0;
348
349 for (int j=0;j<yBasis_.getCardinality();j++) {
350 for (int i=0;i<xBasis_.getCardinality();i++) {
351 for (int k=0;k<inputPoints.dimension(0);k++) {
352 outputValues(bfcur,k,0) = xBasisValues(i,k) * yBasisDerivs(j,k,0);
353 outputValues(bfcur,k,1) = -xBasisDerivs(i,k,0) * yBasisValues(j,k);
354 }
355 bfcur++;
356 }
357 }
358 }
359 break;
360 default:
361 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
362 ">>> ERROR (Basis_HGRAD_QUAD_Cn_FEM): Operator type not implemented");
363 break;
364 }
365 }
366
367 template<class Scalar,class ArrayScalar>
369 const ArrayScalar & inputPoints,
370 const ArrayScalar & cellVertices,
371 const EOperator operatorType) const {
372 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
373 ">>> ERROR (Basis_HGRAD_QUAD_Cn_FEM): FEM Basis calling an FVD member function");
374 }
375
376 template<class Scalar,class ArrayScalar>
378 {
379 int cur = 0;
380 for (int j=0;j<ptsy_.dimension(0);j++)
381 {
382 for (int i=0;i<ptsx_.dimension(0);i++)
383 {
384 dofCoords(cur,0) = ptsx_(i,0);
385 dofCoords(cur,1) = ptsy_(j,0);
386 cur++;
387 }
388 }
389 }
390
391
392}// namespace Intrepid
393
394#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 getDkMultiplicities(Teuchos::Array< int > &partialMult, const int derivativeEnum, const EOperator operatorType, const int spaceDim)
Returns multiplicities of dx, dy, and dz based on the enumeration of the partial derivative,...
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1,...
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
virtual void getDofCoords(ArrayScalar &dofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on a reference cell; defined for interp...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Quadrilateral cell.
Basis_HGRAD_QUAD_Cn_FEM(const int orderx, const int ordery, const ArrayScalar &pts_x, const ArrayScalar &pts_y)
Constructor.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
bool basisTagsAreSet_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
virtual void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const =0
Evaluation of a FEM basis on a reference cell.
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....
virtual int getCardinality() const
Returns cardinality of the basis.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0.
int dimension(const int whichDim) const
Returns the specified dimension.
static void lineProduct2d(const int dim0, const int entity0, const int dim1, const int entity1, int &resultdim, int &resultentity)