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.
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 lineProduct3d(const int dim0, const int entity0, const int dim1, const int entity1, const int dim2, const int entity2, int &resultdim, int &resultentity)