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 
51 namespace 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>
412  void Basis_HGRAD_HEX_Cn_FEM<Scalar, ArrayScalar>::getDofCoords( ArrayScalar & dofCoords ) const
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
virtual int getCardinality() const
Returns cardinality of the basis.
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1...
EBasis basisType_
Type of the basis.
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.
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
bool basisTagsAreSet_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
virtual void getDofCoords(ArrayScalar &DofCoords) const
implement the DofCoordsInterface interface
EOperator
Enumeration of primitive operators available in Intrepid. Primitive operators act on reconstructed fu...
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package http://trilinos...
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
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.
EPointType
Enumeration of types of point distributions in Intrepid.
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)
virtual void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const =0
Evaluation of a FEM basis on a reference cell.
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, its order and the space dimension. Inverse of the getDkEnumeration() method.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Hexahedron cell.
int basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
int basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom. ...