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>
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>
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>
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>
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>
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>
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
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.
static void lineProduct2d(const int dim0, const int entity0, const int dim1, const int entity1, int &resultdim, int &resultentity)
bool basisTagsAreSet_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
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.
Basis_HGRAD_QUAD_Cn_FEM(const int orderx, const int ordery, const ArrayScalar &pts_x, const ArrayScalar &pts_y)
Constructor.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
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.
EPointType
Enumeration of types of point distributions in Intrepid.
virtual void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const =0
Evaluation of a FEM basis on a reference cell.
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
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.
virtual void getDofCoords(ArrayScalar &dofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on a reference cell; defined for interp...
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. ...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Quadrilateral cell.