Intrepid
Intrepid_HDIV_QUAD_In_FEMDef.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Intrepid Package
5 // Copyright (2007) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38 // Denis Ridzal (dridzal@sandia.gov), or
39 // Kara Peterson (kjpeter@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
49 namespace Intrepid {
50 
51  template<class Scalar, class ArrayScalar>
53  const ArrayScalar & ptsClosed ,
54  const ArrayScalar & ptsOpen):
55  closedBasis_( order , ptsClosed ),
56  openBasis_( order-1 , ptsOpen ),
57  closedPts_( ptsClosed ),
58  openPts_( ptsOpen )
59  {
60  this -> basisDegree_ = order;
61  this -> basisCardinality_ = 2 * closedBasis_.getCardinality() * openBasis_.getCardinality();
62  this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
63  this -> basisType_ = BASIS_FEM_FIAT;
64  this -> basisCoordinates_ = COORDINATES_CARTESIAN;
65  this -> basisTagsAreSet_ = false;
66 
67  Array<Array<RCP<Basis<Scalar,ArrayScalar > > > > bases(2);
68  bases[0].resize(2); bases[1].resize(2);
69  bases[0][0] = rcp( &closedBasis_ , false );
70  bases[0][1] = rcp( &openBasis_ , false );
71  bases[1][0] = rcp( &openBasis_ , false );
72  bases[1][1] = rcp( &closedBasis_ , false );
73  this->setBases( bases );
74 
75  }
76 
77  template<class Scalar, class ArrayScalar>
78  Basis_HDIV_QUAD_In_FEM<Scalar,ArrayScalar>::Basis_HDIV_QUAD_In_FEM( int order , const EPointType &pointType ):
79  closedBasis_( order , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL:POINTTYPE_EQUISPACED ),
80  openBasis_( order-1 , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL_OPEN:POINTTYPE_EQUISPACED ),
81  closedPts_( order+1 , 1 ),
82  openPts_( order , 1 )
83  {
84  this -> basisDegree_ = order;
85  this -> basisCardinality_ = 2 * closedBasis_.getCardinality() * openBasis_.getCardinality();
86  this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
87  this -> basisType_ = BASIS_FEM_FIAT;
88  this -> basisCoordinates_ = COORDINATES_CARTESIAN;
89  this -> basisTagsAreSet_ = false;
90 
91  PointTools::getLattice<Scalar,FieldContainer<Scalar> >( closedPts_ ,
92  shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) ,
93  order ,
94  0 ,
95  pointType==POINTTYPE_SPECTRAL?POINTTYPE_WARPBLEND:POINTTYPE_EQUISPACED );
96 
97  if (pointType == POINTTYPE_SPECTRAL)
98  {
99  PointTools::getGaussPoints<Scalar,FieldContainer<Scalar> >( openPts_ ,
100  order - 1 );
101  }
102  else
103  {
104  PointTools::getLattice<Scalar,FieldContainer<Scalar> >( openPts_ ,
105  shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) ,
106  order - 1,
107  0 ,
108  POINTTYPE_EQUISPACED );
109 
110  }
111 
112  Array<Array<RCP<Basis<Scalar,ArrayScalar > > > > bases(2);
113  bases[0].resize(2); bases[1].resize(2);
114  bases[0][0] = rcp( &closedBasis_ , false );
115  bases[0][1] = rcp( &openBasis_ , false );
116  bases[1][0] = rcp( &openBasis_ , false );
117  bases[1][1] = rcp( &closedBasis_ , false );
118  this->setBases( bases );
119 
120  }
121 
122  template<class Scalar, class ArrayScalar>
124 
125  // Basis-dependent intializations
126  int tagSize = 4; // size of DoF tag
127  int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
128  int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
129  int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
130 
131  std::vector<int> tags( tagSize * this->getCardinality() );
132 
133  const std::vector<std::vector<int> >& closedDofTags = closedBasis_.getAllDofTags();
134  const std::vector<std::vector<int> >& openDofTags = openBasis_.getAllDofTags();
135 
136  std::map<int,std::map<int,int> > total_dof_per_entity;
137  std::map<int,std::map<int,int> > current_dof_per_entity;
138 
139  for (int i=0;i<4;i++) {
140  total_dof_per_entity[0][i] = 0;
141  current_dof_per_entity[0][i] = 0;
142  }
143  for (int i=0;i<4;i++) {
144  total_dof_per_entity[1][i] = 0;
145  current_dof_per_entity[1][i] = 0;
146  }
147  total_dof_per_entity[2][0] = 0;
148  current_dof_per_entity[2][0] = 0;
149 
150  // tally dof on each facet. none on vertex
151  for (int i=0;i<4;i++) {
152  total_dof_per_entity[1][i] = openBasis_.getCardinality();
153  }
154 
155  total_dof_per_entity[2][0] = this->getCardinality() - 4 * openBasis_.getCardinality();
156 
157  int tagcur = 0;
158  // loop over the x-component basis functions, which are (psi(x)phi(y),0)
159  // for psi in the closed basis and phi in the open
160  for (int j=0;j<openBasis_.getCardinality();j++) {
161  const int odim = openDofTags[j][0];
162  const int oent = openDofTags[j][1];
163  for (int i=0;i<closedBasis_.getCardinality();i++) {
164  const int cdim = closedDofTags[i][0];
165  const int cent = closedDofTags[i][1];
166  int dofdim;
167  int dofent;
168  ProductTopology::lineProduct2d(cdim,cent,odim,oent,dofdim,dofent);
169  tags[4*tagcur] = dofdim;
170  tags[4*tagcur+1] = dofent;
171  tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
172  current_dof_per_entity[dofdim][dofent]++;
173  tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
174  tagcur++;
175  }
176  }
177  // now we have to do it for the y-component basis functions, which are
178  // (0,phi(x)psi(y)) for psi in the closed basis and phi in the open
179  for (int j=0;j<closedBasis_.getCardinality();j++) {
180  const int cdim = closedDofTags[j][0];
181  const int cent = closedDofTags[j][1];
182  for (int i=0;i<openBasis_.getCardinality();i++) {
183  const int odim = openDofTags[i][0];
184  const int oent = openDofTags[i][1];
185  int dofdim;
186  int dofent;
187  ProductTopology::lineProduct2d(odim,oent,cdim,cent,dofdim,dofent);
188  tags[4*tagcur] = dofdim;
189  tags[4*tagcur+1] = dofent;
190  tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
191  current_dof_per_entity[dofdim][dofent]++;
192  tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
193  tagcur++;
194  }
195  }
196 
197 // for (int i=0;i<this->getCardinality();i++) {
198 // for (int j=0;j<4;j++) {
199 // std::cout << tags[4*i+j] << " ";
200 // }
201 // std::cout << std::endl;
202 // }
203 
204  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
205  Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
206  this -> ordinalToTag_,
207  &(tags[0]),
208  this -> basisCardinality_,
209  tagSize,
210  posScDim,
211  posScOrd,
212  posDfOrd);
213  }
214 
215 
216  template<class Scalar, class ArrayScalar>
218  const ArrayScalar & inputPoints,
219  const EOperator operatorType) const {
220 
221  // Verify arguments
222 #ifdef HAVE_INTREPID_DEBUG
223  Intrepid::getValues_HDIV_Args<Scalar, ArrayScalar>(outputValues,
224  inputPoints,
225  operatorType,
226  this -> getBaseCellTopology(),
227  this -> getCardinality() );
228 #endif
229 
230  // Number of evaluation points = dim 0 of inputPoints
231  int dim0 = inputPoints.dimension(0);
232 
233  // separate out points
234  FieldContainer<Scalar> xPoints(dim0,1);
235  FieldContainer<Scalar> yPoints(dim0,1);
236 
237  for (int i=0;i<dim0;i++) {
238  xPoints(i,0) = inputPoints(i,0);
239  yPoints(i,0) = inputPoints(i,1);
240  }
241 
242  switch (operatorType) {
243  case OPERATOR_VALUE:
244  {
245  FieldContainer<Scalar> closedBasisValsXPts( closedBasis_.getCardinality() , dim0 );
246  FieldContainer<Scalar> closedBasisValsYPts( closedBasis_.getCardinality() , dim0 );
247  FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 );
248  FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 );
249 
250  closedBasis_.getValues( closedBasisValsXPts , xPoints , OPERATOR_VALUE );
251  closedBasis_.getValues( closedBasisValsYPts , yPoints , OPERATOR_VALUE );
252  openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE );
253  openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE );
254 
255  int bfcur = 0;
256  // x component bfs are (closed(x) open(y),0)
257  for (int j=0;j<openBasis_.getCardinality();j++) {
258  for (int i=0;i<closedBasis_.getCardinality();i++) {
259  for (int l=0;l<dim0;l++) {
260  outputValues(bfcur,l,0) = closedBasisValsXPts(i,l) * openBasisValsYPts(j,l);
261  outputValues(bfcur,l,1) = 0.0;
262  }
263  bfcur++;
264  }
265  }
266 
267  // y component bfs are (0,open(x) closed(y))
268  for (int j=0;j<closedBasis_.getCardinality();j++) {
269  for (int i=0;i<openBasis_.getCardinality();i++) {
270  for (int l=0;l<dim0;l++) {
271  outputValues(bfcur,l,0) = 0.0;
272  outputValues(bfcur,l,1) = openBasisValsXPts(i,l) * closedBasisValsYPts(j,l);
273  }
274  bfcur++;
275  }
276  }
277  }
278  break;
279  case OPERATOR_DIV:
280  {
281  FieldContainer<Scalar> closedBasisDerivsXPts( closedBasis_.getCardinality() , dim0 , 1 );
282  FieldContainer<Scalar> closedBasisDerivsYPts( closedBasis_.getCardinality() , dim0 , 1 );
283  FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 );
284  FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 );
285 
286  closedBasis_.getValues( closedBasisDerivsXPts , xPoints , OPERATOR_D1 );
287  closedBasis_.getValues( closedBasisDerivsYPts , yPoints , OPERATOR_D1 );
288  openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE );
289  openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE );
290 
291  int bfcur = 0;
292 
293  // x component basis functions first
294  for (int j=0;j<openBasis_.getCardinality();j++) {
295  for (int i=0;i<closedBasis_.getCardinality();i++) {
296  for (int l=0;l<dim0;l++) {
297  outputValues(bfcur,l) = closedBasisDerivsXPts(i,l,0) * openBasisValsYPts(j,l);
298  }
299  bfcur++;
300  }
301  }
302 
303  // now y component basis functions
304  for (int j=0;j<closedBasis_.getCardinality();j++) {
305  for (int i=0;i<openBasis_.getCardinality();i++) {
306  for (int l=0;l<dim0;l++) {
307  outputValues(bfcur,l) = openBasisValsXPts(i,l) * closedBasisDerivsYPts(j,l,0);
308  }
309  bfcur++;
310  }
311  }
312  }
313  break;
314  case OPERATOR_CURL:
315  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument,
316  ">>> ERROR (Basis_HDIV_QUAD_In_FEM): CURL is invalid operator for HDIV Basis Functions");
317  break;
318 
319  case OPERATOR_GRAD:
320  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_GRAD), std::invalid_argument,
321  ">>> ERROR (Basis_HDIV_QUAD_In_FEM): GRAD is invalid operator for HDIV Basis Functions");
322  break;
323 
324  case OPERATOR_D1:
325  case OPERATOR_D2:
326  case OPERATOR_D3:
327  case OPERATOR_D4:
328  case OPERATOR_D5:
329  case OPERATOR_D6:
330  case OPERATOR_D7:
331  case OPERATOR_D8:
332  case OPERATOR_D9:
333  case OPERATOR_D10:
334  TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType == OPERATOR_D1) ||
335  (operatorType == OPERATOR_D2) ||
336  (operatorType == OPERATOR_D3) ||
337  (operatorType == OPERATOR_D4) ||
338  (operatorType == OPERATOR_D5) ||
339  (operatorType == OPERATOR_D6) ||
340  (operatorType == OPERATOR_D7) ||
341  (operatorType == OPERATOR_D8) ||
342  (operatorType == OPERATOR_D9) ||
343  (operatorType == OPERATOR_D10) ),
344  std::invalid_argument,
345  ">>> ERROR (Basis_HDIV_QUAD_In_FEM): Invalid operator type");
346  break;
347 
348  default:
349  TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType != OPERATOR_VALUE) &&
350  (operatorType != OPERATOR_GRAD) &&
351  (operatorType != OPERATOR_CURL) &&
352  (operatorType != OPERATOR_DIV) &&
353  (operatorType != OPERATOR_D1) &&
354  (operatorType != OPERATOR_D2) &&
355  (operatorType != OPERATOR_D3) &&
356  (operatorType != OPERATOR_D4) &&
357  (operatorType != OPERATOR_D5) &&
358  (operatorType != OPERATOR_D6) &&
359  (operatorType != OPERATOR_D7) &&
360  (operatorType != OPERATOR_D8) &&
361  (operatorType != OPERATOR_D9) &&
362  (operatorType != OPERATOR_D10) ),
363  std::invalid_argument,
364  ">>> ERROR (Basis_HDIV_QUAD_In_FEM): Invalid operator type");
365  }
366  }
367 
368 
369 
370  template<class Scalar, class ArrayScalar>
372  const ArrayScalar & inputPoints,
373  const ArrayScalar & cellVertices,
374  const EOperator operatorType) const {
375  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
376  ">>> ERROR (Basis_HDIV_QUAD_In_FEM): FEM Basis calling an FVD member function");
377  }
378 
379  template<class Scalar, class ArrayScalar>
381  {
382  // x-component basis functions
383  int cur = 0;
384 
385  for (int j=0;j<openPts_.dimension(0);j++)
386  {
387  for (int i=0;i<closedPts_.dimension(0);i++)
388  {
389  DofCoords(cur,0) = closedPts_(i,0);
390  DofCoords(cur,1) = openPts_(j,0);
391  cur++;
392  }
393  }
394 
395  // y-component basis functions
396  for (int j=0;j<closedPts_.dimension(0);j++)
397  {
398  for (int i=0;i<openPts_.dimension(0);i++)
399  {
400  DofCoords(cur,0) = openPts_(i,0);
401  DofCoords(cur,1) = closedPts_(j,0);
402  cur++;
403  }
404  }
405 
406  return;
407  }
408 
409 
410 }// namespace Intrepid
virtual void getDofCoords(ArrayScalar &DofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on a reference cell; defined for interp...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Quadrilateral cell.
EBasis basisType_
Type of the basis.
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
static void lineProduct2d(const int dim0, const int entity0, const int dim1, const int entity1, int &resultdim, int &resultentity)
bool basisTagsAreSet_
&quot;true&quot; if tagToOrdinal_ and ordinalToTag_ have been initialized
Basis_HDIV_QUAD_In_FEM(int order, const ArrayScalar &ptsClosed, const ArrayScalar &ptsOpen)
Constructor.
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package http://trilinos...
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. ...