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 
51 using Teuchos::Array;
52 using Teuchos::RCP;
53 
54 namespace Intrepid {
55  template<class Scalar, class ArrayScalar>
56  Basis_HGRAD_QUAD_Cn_FEM<Scalar,ArrayScalar>::Basis_HGRAD_QUAD_Cn_FEM( const int orderx , const int ordery,
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>
93  Basis_HGRAD_QUAD_Cn_FEM<Scalar,ArrayScalar>::Basis_HGRAD_QUAD_Cn_FEM( const int order,
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>
126  void Basis_HGRAD_QUAD_Cn_FEM<Scalar,ArrayScalar>::initializeTags()
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>
213  void Basis_HGRAD_QUAD_Cn_FEM<Scalar,ArrayScalar>::getValues( ArrayScalar &outputValues ,
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>
368  void Basis_HGRAD_QUAD_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues,
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>
377  void Basis_HGRAD_QUAD_Cn_FEM<Scalar, ArrayScalar>::getDofCoords( ArrayScalar & dofCoords ) const
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
395 
396 #if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
397 #ifdef __GNUC__
398 #warning "The Intrepid package is deprecated"
399 #endif
400 #endif
401 
Intrepid::Basis::getCardinality
virtual int getCardinality() const
Returns cardinality of the basis.
Definition: Intrepid_BasisDef.hpp:100
Intrepid::Basis_HGRAD_LINE_Cn_FEM
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1...
Definition: Intrepid_HGRAD_LINE_Cn_FEM.hpp:78
Intrepid::Basis::basisType_
EBasis basisType_
Type of the basis.
Definition: Intrepid_Basis.hpp:113
Intrepid::Basis::getAllDofTags
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
Definition: Intrepid_BasisDef.hpp:89
Intrepid::ProductTopology::lineProduct2d
static void lineProduct2d(const int dim0, const int entity0, const int dim1, const int entity1, int &resultdim, int &resultentity)
Definition: Intrepid_ProductTopology.hpp:75
Intrepid::Basis::basisTagsAreSet_
bool basisTagsAreSet_
&quot;true&quot; if tagToOrdinal_ and ordinalToTag_ have been initialized
Definition: Intrepid_Basis.hpp:121
Intrepid::Basis::basisCoordinates_
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
Definition: Intrepid_Basis.hpp:117
Intrepid::Basis_HGRAD_QUAD_Cn_FEM::Basis_HGRAD_QUAD_Cn_FEM
Basis_HGRAD_QUAD_Cn_FEM(const int orderx, const int ordery, const ArrayScalar &pts_x, const ArrayScalar &pts_y)
Constructor.
Definition: Intrepid_HGRAD_QUAD_Cn_FEMDef.hpp:56
Intrepid::Basis
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
Definition: Intrepid_Basis.hpp:89
Intrepid::Basis_HGRAD_QUAD_Cn_FEM::getValues
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Quadrilateral cell.
Definition: Intrepid_HGRAD_QUAD_Cn_FEMDef.hpp:213
Intrepid::Basis::basisCellTopology_
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package http://trilinos...
Definition: Intrepid_Basis.hpp:109
Intrepid::Basis_HGRAD_QUAD_Cn_FEM::getDofCoords
virtual void getDofCoords(ArrayScalar &dofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on a reference cell; defined for interp...
Definition: Intrepid_HGRAD_QUAD_Cn_FEMDef.hpp:377
Intrepid::Basis::getValues
virtual void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const =0
Evaluation of a FEM basis on a reference cell.
Intrepid::Basis_HGRAD_QUAD_Cn_FEM::initializeTags
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
Definition: Intrepid_HGRAD_QUAD_Cn_FEMDef.hpp:126
Intrepid::Basis::basisDegree_
int basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
Definition: Intrepid_Basis.hpp:104
Intrepid::Basis::basisCardinality_
int basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom. ...
Definition: Intrepid_Basis.hpp:100
Generated by   doxygen 1.8.5