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...
virtual void getDofCoords(ArrayScalar &DofCoords) const
implement the DofCoordsInterface interface
EBasis basisType_
Type of the basis.
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
bool basisTagsAreSet_
&quot;true&quot; if tagToOrdinal_ and ordinalToTag_ have been initialized
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...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Hexahedron cell.
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.
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.
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. ...