Intrepid
Intrepid_HGRAD_HEX_Cn_FEMDef.hpp
Go to the documentation of this file.
1#ifndef INTREPID_HGRAD_HEX_CN_FEMDEF_HPP
2#define INTREPID_HGRAD_HEX_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 template<class Scalar, class ArrayScalar>
54 const int ordery ,
55 const int orderz ,
56 const ArrayScalar &pts_x ,
57 const ArrayScalar &pts_y ,
58 const ArrayScalar &pts_z ):
59 ptsx_( pts_x.dimension(0) , 1 ),
60 ptsy_( pts_y.dimension(0) , 1 ),
61 ptsz_( pts_z.dimension(0) , 1 )
62 {
63 for (int i=0;i<pts_x.dimension(0);i++)
64 {
65 ptsx_(i,0) = pts_x(i,0);
66 }
67 for (int i=0;i<pts_y.dimension(0);i++)
68 {
69 ptsy_(i,0) = pts_y(i,0);
70 }
71 for (int i=0;i<pts_z.dimension(0);i++)
72 {
73 ptsz_(i,0) = pts_z(i,0);
74 }
75
76 Array<Array<RCP<Basis<Scalar,ArrayScalar> > > > bases(1);
77 bases[0].resize(3);
78
79 bases[0][0] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM< Scalar , ArrayScalar >( orderx , pts_x ) );
80 bases[0][1] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM< Scalar , ArrayScalar >( ordery , pts_y ) );
81 bases[0][2] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM< Scalar , ArrayScalar >( orderz , pts_z ) );
82
83 this->setBases( bases );
84
85 this->basisCardinality_ = (orderx+1)*(ordery+1)*(orderz+1);
86 if (orderx >= ordery && orderx >= orderz ) {
87 this->basisDegree_ = orderx;
88 }
89 else if (ordery >= orderx && ordery >= orderz) {
90 this->basisDegree_ = ordery;
91 }
92 else {
93 this->basisDegree_ = orderz;
94 }
95 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() );
96 this -> basisType_ = BASIS_FEM_FIAT;
97 this -> basisCoordinates_ = COORDINATES_CARTESIAN;
98 this -> basisTagsAreSet_ = false;
99 }
100
101 template<class Scalar, class ArrayScalar>
103 const EPointType & pointType ):
104 ptsx_( order+1 , 1 ),
105 ptsy_( order+1 , 1 ),
106 ptsz_( order+1 , 1 )
107 {
108 Array<Array<RCP<Basis<Scalar,ArrayScalar> > > > bases(1);
109 bases[0].resize(3);
110
111 bases[0][0] = Teuchos::rcp( new Basis_HGRAD_LINE_Cn_FEM< Scalar , ArrayScalar >( order , pointType ) );
112 // basis is same in each direction, so I only need to instantiate it once!
113 bases[0][1] = bases[0][0];
114 bases[0][2] = bases[0][0];
115
116 this->setBases( bases );
117
118 this->basisCardinality_ = (order+1)*(order+1)*(order+1);
119 this->basisDegree_ = order;
120 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() );
121 this -> basisType_ = BASIS_FEM_FIAT;
122 this -> basisCoordinates_ = COORDINATES_CARTESIAN;
123 this -> basisTagsAreSet_ = false;
124
125 // get points
126 EPointType pt = (pointType==POINTTYPE_EQUISPACED)?pointType:POINTTYPE_WARPBLEND;
127 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( ptsx_ ,
128 shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) ,
129 order ,
130 0 ,
131 pt );
132 for (int i=0;i<order+1;i++)
133 {
134 ptsy_(i,0) = ptsx_(i,0);
135 ptsz_(i,0) = ptsx_(i,0);
136 }
137
138 }
139
140 template<class Scalar, class ArrayScalar>
142 {
143 // Basis-dependent initializations
144 int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
145 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
146 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
147 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
148
149 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
150
151 std::vector<int> tags( tagSize * this->getCardinality() );
152
153 // temporarily just put everything on the cell itself
154 for (int i=0;i<this->getCardinality();i++) {
155 tags[tagSize*i] = 3;
156 tags[tagSize*i+1] = 0;
157 tags[tagSize*i+2] = i;
158 tags[tagSize*i+3] = this->getCardinality();
159 }
160
161 Basis<Scalar,ArrayScalar> &xBasis_ = *this->bases_[0][0];
162 Basis<Scalar,ArrayScalar> &yBasis_ = *this->bases_[0][1];
163 Basis<Scalar,ArrayScalar> &zBasis_ = *this->bases_[0][2];
164
165
166 // now let's try to do it "right"
167 // let's get the x, y, z bases and their dof
168 const std::vector<std::vector<int> >& xdoftags = xBasis_.getAllDofTags();
169 const std::vector<std::vector<int> >& ydoftags = yBasis_.getAllDofTags();
170 const std::vector<std::vector<int> >& zdoftags = zBasis_.getAllDofTags();
171
172 std::map<int,std::map<int,int> > total_dof_per_entity;
173 std::map<int,std::map<int,int> > current_dof_per_entity;
174
175 // vertex dof
176 for (int i=0;i<8;i++) {
177 total_dof_per_entity[0][i] = 0;
178 current_dof_per_entity[0][i] = 0;
179 }
180 // edge dof (12 edges)
181 for (int i=0;i<12;i++) {
182 total_dof_per_entity[1][i] = 0;
183 current_dof_per_entity[1][i] = 0;
184 }
185 // face dof (6 faces)
186 for (int i=0;i<6;i++) {
187 total_dof_per_entity[2][i] = 0;
188 current_dof_per_entity[2][i] = 0;
189 }
190 // internal dof
191 total_dof_per_entity[3][0] = 0;
192 //total_dof_per_entity[3][1] = 0;
193
194 // let's tally the total degrees of freedom on each entity
195 for (int k=0;k<zBasis_.getCardinality();k++) {
196 const int zdim = zdoftags[k][0];
197 const int zent = zdoftags[k][1];
198 for (int j=0;j<yBasis_.getCardinality();j++) {
199 const int ydim = ydoftags[j][0];
200 const int yent = ydoftags[j][1];
201 for (int i=0;i<xBasis_.getCardinality();i++) {
202 const int xdim = xdoftags[i][0];
203 const int xent = xdoftags[i][1];
204 int dofdim;
205 int dofent;
206 ProductTopology::lineProduct3d( xdim , xent , ydim , yent , zdim , zent , dofdim , dofent );
207 total_dof_per_entity[dofdim][dofent] += 1;
208
209 }
210 }
211 }
212
213 int tagcur = 0;
214 for (int k=0;k<zBasis_.getCardinality();k++) {
215 const int zdim = zdoftags[k][0];
216 const int zent = zdoftags[k][1];
217 for (int j=0;j<yBasis_.getCardinality();j++) {
218 const int ydim = ydoftags[j][0];
219 const int yent = ydoftags[j][1];
220 for (int i=0;i<xBasis_.getCardinality();i++) {
221 const int xdim = xdoftags[i][0];
222 const int xent = xdoftags[i][1];
223 int dofdim;
224 int dofent;
225 ProductTopology::lineProduct3d( xdim , xent , ydim , yent , zdim , zent , dofdim , dofent );
226 tags[4*tagcur] = dofdim;
227 tags[4*tagcur+1] = dofent;
228 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
229 current_dof_per_entity[dofdim][dofent]++;
230 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
231 tagcur++;
232 }
233 }
234 }
235
236 Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
237 this -> ordinalToTag_,
238 &(tags[0]),
239 this -> basisCardinality_,
240 tagSize,
241 posScDim,
242 posScOrd,
243 posDfOrd);
244
245 }
246
247 template<class Scalar, class ArrayScalar>
249 const ArrayScalar &inputPoints ,
250 const EOperator operatorType ) const
251 {
252#ifdef HAVE_INTREPID_DEBUG
253 Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
254 inputPoints,
255 operatorType,
256 this -> getBaseCellTopology(),
257 this -> getCardinality() );
258#endif
259
260 Basis<Scalar,ArrayScalar> &xBasis_ = *this->bases_[0][0];
261 Basis<Scalar,ArrayScalar> &yBasis_ = *this->bases_[0][1];
262 Basis<Scalar,ArrayScalar> &zBasis_ = *this->bases_[0][2];
263
264
265 FieldContainer<Scalar> xInputPoints(inputPoints.dimension(0),1);
266 FieldContainer<Scalar> yInputPoints(inputPoints.dimension(0),1);
267 FieldContainer<Scalar> zInputPoints(inputPoints.dimension(0),1);
268
269 for (int i=0;i<inputPoints.dimension(0);i++) {
270 xInputPoints(i,0) = inputPoints(i,0);
271 yInputPoints(i,0) = inputPoints(i,1);
272 zInputPoints(i,0) = inputPoints(i,2);
273 }
274
275 switch (operatorType) {
276 case OPERATOR_VALUE:
277 {
278 FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0));
279 FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0));
280 FieldContainer<Scalar> zBasisValues(zBasis_.getCardinality(),zInputPoints.dimension(0));
281
282 xBasis_.getValues(xBasisValues,xInputPoints,OPERATOR_VALUE);
283 yBasis_.getValues(yBasisValues,yInputPoints,OPERATOR_VALUE);
284 zBasis_.getValues(zBasisValues,zInputPoints,OPERATOR_VALUE);
285
286 int bfcur = 0;
287 for (int k=0;k<zBasis_.getCardinality();k++) {
288 for (int j=0;j<yBasis_.getCardinality();j++) {
289 for (int i=0;i<xBasis_.getCardinality();i++) {
290 for (int l=0;l<inputPoints.dimension(0);l++) {
291 outputValues(bfcur,l) = xBasisValues(i,l) * yBasisValues(j,l) * zBasisValues(k,l);
292 }
293 bfcur++;
294 }
295 }
296 }
297 }
298 break;
299 case OPERATOR_D1:
300 case OPERATOR_GRAD:
301 {
302 FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0));
303 FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0));
304 FieldContainer<Scalar> zBasisValues(zBasis_.getCardinality(),zInputPoints.dimension(0));
305 FieldContainer<Scalar> xBasisDerivs(xBasis_.getCardinality(),xInputPoints.dimension(0),1);
306 FieldContainer<Scalar> yBasisDerivs(yBasis_.getCardinality(),yInputPoints.dimension(0),1);
307 FieldContainer<Scalar> zBasisDerivs(zBasis_.getCardinality(),zInputPoints.dimension(0),1);
308
309 xBasis_.getValues(xBasisValues,xInputPoints,OPERATOR_VALUE);
310 yBasis_.getValues(yBasisValues,yInputPoints,OPERATOR_VALUE);
311 zBasis_.getValues(zBasisValues,zInputPoints,OPERATOR_VALUE);
312 xBasis_.getValues(xBasisDerivs,xInputPoints,OPERATOR_D1);
313 yBasis_.getValues(yBasisDerivs,yInputPoints,OPERATOR_D1);
314 zBasis_.getValues(zBasisDerivs,zInputPoints,OPERATOR_D1);
315
316 int bfcur = 0;
317 for (int k=0;k<zBasis_.getCardinality();k++) {
318 for (int j=0;j<yBasis_.getCardinality();j++) {
319 for (int i=0;i<xBasis_.getCardinality();i++) {
320 for (int l=0;l<inputPoints.dimension(0);l++) {
321 outputValues(bfcur,l,0) = xBasisDerivs(i,l,0) * yBasisValues(j,l) * zBasisValues(k,l);
322 outputValues(bfcur,l,1) = xBasisValues(i,l) * yBasisDerivs(j,l,0) * zBasisValues(k,l);
323 outputValues(bfcur,l,2) = xBasisValues(i,l) * yBasisValues(j,l) * zBasisDerivs(k,l,0);
324 }
325 bfcur++;
326 }
327 }
328 }
329 }
330 break;
331 case OPERATOR_D2:
332 case OPERATOR_D3:
333 case OPERATOR_D4:
334 case OPERATOR_D5:
335 case OPERATOR_D6:
336 case OPERATOR_D7:
337 case OPERATOR_D8:
338 case OPERATOR_D9:
339 case OPERATOR_D10:
340 {
341 FieldContainer<Scalar> xBasisValues(xBasis_.getCardinality(),xInputPoints.dimension(0));
342 FieldContainer<Scalar> yBasisValues(yBasis_.getCardinality(),yInputPoints.dimension(0));
343 FieldContainer<Scalar> zBasisValues(yBasis_.getCardinality(),zInputPoints.dimension(0));
344
345 Teuchos::Array<int> partialMult;
346
347 for (int d=0;d<getDkCardinality(operatorType,3);d++) {
348 getDkMultiplicities( partialMult , d , operatorType , 3 );
349 if (partialMult[0] == 0) {
350 xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0));
351 xBasis_.getValues( xBasisValues , xInputPoints, OPERATOR_VALUE );
352 }
353 else {
354 xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0),1);
355 EOperator xop = (EOperator) ( (int) OPERATOR_D1 + partialMult[0] - 1 );
356 xBasis_.getValues( xBasisValues , xInputPoints, xop );
357 xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0));
358 }
359 if (partialMult[1] == 0) {
360 yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0));
361 yBasis_.getValues( yBasisValues , yInputPoints, OPERATOR_VALUE );
362 }
363 else {
364 yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0),1);
365 EOperator yop = (EOperator) ( (int) OPERATOR_D1 + partialMult[1] - 1 );
366 yBasis_.getValues( yBasisValues , yInputPoints, yop );
367 yBasisValues.resize(yBasis_.getCardinality(),yInputPoints.dimension(0));
368 }
369 if (partialMult[2] == 0) {
370 zBasisValues.resize(zBasis_.getCardinality(),zInputPoints.dimension(0));
371 zBasis_.getValues( zBasisValues , zInputPoints, OPERATOR_VALUE );
372 }
373 else {
374 zBasisValues.resize(zBasis_.getCardinality(),zInputPoints.dimension(0),1);
375 EOperator zop = (EOperator) ( (int) OPERATOR_D1 + partialMult[2] - 1 );
376 zBasis_.getValues( zBasisValues , zInputPoints, zop );
377 zBasisValues.resize(zBasis_.getCardinality(),zInputPoints.dimension(0));
378 }
379
380
381 int bfcur = 0;
382 for (int k=0;k<zBasis_.getCardinality();k++) {
383 for (int j=0;j<yBasis_.getCardinality();j++) {
384 for (int i=0;i<xBasis_.getCardinality();i++) {
385 for (int l=0;l<inputPoints.dimension(0);l++) {
386 outputValues(bfcur,l,d) = xBasisValues(i,l) * yBasisValues(j,l) * zBasisValues(k,l);
387 }
388 bfcur++;
389 }
390 }
391 }
392 }
393 }
394 break;
395 default:
396 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
397 ">>> ERROR (Basis_HGRAD_HEX_Cn_FEM): Operator type not implemented");
398 break;
399 }
400 }
401
402 template<class Scalar,class ArrayScalar>
404 const ArrayScalar & inputPoints,
405 const ArrayScalar & cellVertices,
406 const EOperator operatorType) const {
407 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
408 ">>> ERROR (Basis_HGRAD_HEX_Cn_FEM): FEM Basis calling an FVD member function");
409}
410
411 template<class Scalar,class ArrayScalar>
413 {
414 int cur = 0;
415 for (int k=0;k<ptsz_.dimension(0);k++)
416 {
417 for (int j=0;j<ptsy_.dimension(0);j++)
418 {
419 for (int i=0;i<ptsx_.dimension(0);i++)
420 {
421 dofCoords(cur,0) = ptsx_(i,0);
422 dofCoords(cur,1) = ptsy_(j,0);
423 dofCoords(cur,2) = ptsz_(k,0);
424 cur++;
425 }
426 }
427 }
428 }
429
430}// namespace Intrepid
431
432#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,...
Basis_HGRAD_HEX_Cn_FEM(const int orderx, const int ordery, const int orderz, const ArrayScalar &pts_x, const ArrayScalar &pts_y, const ArrayScalar &pts_z)
Constructor.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Hexahedron cell.
virtual void getDofCoords(ArrayScalar &DofCoords) const
implement the DofCoordsInterface interface
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1,...
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 lineProduct3d(const int dim0, const int entity0, const int dim1, const int entity1, const int dim2, const int entity2, int &resultdim, int &resultentity)