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,...
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.
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....
static void lineProduct2d(const int dim0, const int entity0, const int dim1, const int entity1, int &resultdim, int &resultentity)