Intrepid
Intrepid_HCURL_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( &openBasis_ , false );
70  bases[0][1] = rcp( &closedBasis_ , false );
71  bases[1][0] = rcp( &closedBasis_ , false );
72  bases[1][1] = rcp( &openBasis_ , false );
73  this->setBases( bases );
74 
75 
76  }
77 
78  template<class Scalar, class ArrayScalar>
80  closedBasis_( order , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL:POINTTYPE_EQUISPACED ),
81  openBasis_( order-1 , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL_OPEN:POINTTYPE_EQUISPACED ) ,
82  closedPts_( order+1 , 1 ),
83  openPts_( order , 1 )
84  {
85  this -> basisDegree_ = order;
86  this -> basisCardinality_ = 2 * closedBasis_.getCardinality() * openBasis_.getCardinality();
87  this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
88  this -> basisType_ = BASIS_FEM_FIAT;
89  this -> basisCoordinates_ = COORDINATES_CARTESIAN;
90  this -> basisTagsAreSet_ = false;
91 
92  PointTools::getLattice<Scalar,FieldContainer<Scalar> >( closedPts_ ,
93  shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) ,
94  order ,
95  0 ,
96  pointType==POINTTYPE_SPECTRAL?POINTTYPE_WARPBLEND:POINTTYPE_EQUISPACED );
97 
98  if (pointType == POINTTYPE_SPECTRAL)
99  {
100  PointTools::getGaussPoints<Scalar,FieldContainer<Scalar> >( openPts_ ,
101  order - 1 );
102  }
103  else
104  {
105  PointTools::getLattice<Scalar,FieldContainer<Scalar> >( openPts_ ,
106  shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) ,
107  order + 1,
108  1 ,
109  POINTTYPE_EQUISPACED );
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( &openBasis_ , false );
115  bases[0][1] = rcp( &closedBasis_ , false );
116  bases[1][0] = rcp( &closedBasis_ , false );
117  bases[1][1] = rcp( &openBasis_ , 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-face tangent basis functions, which are (phi(x)psi(y),0)
159  // for psi in the closed basis and phi in the open
160  for (int j=0;j<closedBasis_.getCardinality();j++)
161  {
162  const int cdim = closedDofTags[j][0];
163  const int cent = closedDofTags[j][1];
164  for (int i=0;i<openBasis_.getCardinality();i++)
165  {
166  const int odim = openDofTags[i][0];
167  const int oent = openDofTags[i][1];
168  int dofdim;
169  int dofent;
170  ProductTopology::lineProduct2d(odim,oent,cdim,cent,dofdim,dofent);
171  tags[4*tagcur] = dofdim;
172  tags[4*tagcur+1] = dofent;
173  tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
174  current_dof_per_entity[dofdim][dofent]++;
175  tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
176  tagcur++;
177  }
178  }
179 
180  // now we have to do it for the other basis functions, which are
181  // (psi(x)phi(y)) for psi in the closed basis and phi in the open
182  for (int j=0;j<openBasis_.getCardinality();j++)
183  {
184  const int odim = openDofTags[j][0];
185  const int oent = openDofTags[j][1];
186  for (int i=0;i<closedBasis_.getCardinality();i++)
187  {
188  const int cdim = closedDofTags[i][0];
189  const int cent = closedDofTags[i][1];
190  int dofdim;
191  int dofent;
192  ProductTopology::lineProduct2d(cdim,cent,odim,oent,dofdim,dofent);
193  tags[4*tagcur] = dofdim;
194  tags[4*tagcur+1] = dofent;
195  tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
196  current_dof_per_entity[dofdim][dofent]++;
197  tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
198  tagcur++;
199  }
200  }
201 
202  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
203  Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
204  this -> ordinalToTag_,
205  &(tags[0]),
206  this -> basisCardinality_,
207  tagSize,
208  posScDim,
209  posScOrd,
210  posDfOrd);
211  }
212 
213 
214  template<class Scalar, class ArrayScalar>
216  const ArrayScalar & inputPoints,
217  const EOperator operatorType) const {
218 
219  // Verify arguments
220 #ifdef HAVE_INTREPID_DEBUG
221  Intrepid::getValues_HCURL_Args<Scalar, ArrayScalar>(outputValues,
222  inputPoints,
223  operatorType,
224  this -> getBaseCellTopology(),
225  this -> getCardinality() );
226 #endif
227 
228  // Number of evaluation points = dim 0 of inputPoints
229  int dim0 = inputPoints.dimension(0);
230 
231  // separate out points
232  FieldContainer<Scalar> xPoints(dim0,1);
233  FieldContainer<Scalar> yPoints(dim0,1);
234 
235  for (int i=0;i<dim0;i++) {
236  xPoints(i,0) = inputPoints(i,0);
237  yPoints(i,0) = inputPoints(i,1);
238  }
239 
240  switch (operatorType) {
241  case OPERATOR_VALUE:
242  {
243  FieldContainer<Scalar> closedBasisValsXPts( closedBasis_.getCardinality() , dim0 );
244  FieldContainer<Scalar> closedBasisValsYPts( closedBasis_.getCardinality() , dim0 );
245  FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 );
246  FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 );
247 
248  closedBasis_.getValues( closedBasisValsXPts , xPoints , OPERATOR_VALUE );
249  closedBasis_.getValues( closedBasisValsYPts , yPoints , OPERATOR_VALUE );
250  openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE );
251  openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE );
252 
253  int bfcur = 0;
254  // x component bfs are (open(x) closed(y),0)
255  for (int j=0;j<closedBasis_.getCardinality();j++)
256  {
257  for (int i=0;i<openBasis_.getCardinality();i++)
258  {
259  for (int l=0;l<dim0;l++)
260  {
261  outputValues(bfcur,l,0) = closedBasisValsYPts(j,l) * openBasisValsXPts(i,l);
262  outputValues(bfcur,l,1) = 0.0;
263  }
264  bfcur++;
265  }
266  }
267 
268  // y component bfs are (0,closed(x) open(y))
269  for (int j=0;j<openBasis_.getCardinality();j++)
270  {
271  for (int i=0;i<closedBasis_.getCardinality();i++)
272  {
273  for (int l=0;l<dim0;l++)
274  {
275  outputValues(bfcur,l,0) = 0.0;
276  outputValues(bfcur,l,1) = openBasisValsYPts(j,l) * closedBasisValsXPts(i,l);
277  }
278  bfcur++;
279  }
280  }
281  }
282  break;
283  case OPERATOR_CURL:
284  {
285  FieldContainer<Scalar> closedBasisDerivsXPts( closedBasis_.getCardinality() , dim0 , 1 );
286  FieldContainer<Scalar> closedBasisDerivsYPts( closedBasis_.getCardinality() , dim0 , 1 );
287  FieldContainer<Scalar> openBasisValsXPts( openBasis_.getCardinality() , dim0 );
288  FieldContainer<Scalar> openBasisValsYPts( openBasis_.getCardinality() , dim0 );
289 
290  closedBasis_.getValues( closedBasisDerivsXPts , xPoints , OPERATOR_D1 );
291  closedBasis_.getValues( closedBasisDerivsYPts , yPoints , OPERATOR_D1 );
292  openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE );
293  openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE );
294 
295  int bfcur = 0;
296 
297  // x component basis functions first
298  for (int j=0;j<closedBasis_.getCardinality();j++)
299  {
300  for (int i=0;i<openBasis_.getCardinality();i++)
301  {
302  for (int l=0;l<dim0;l++) {
303  outputValues(bfcur,l) = -closedBasisDerivsYPts(j,l,0) * openBasisValsXPts(i,l);
304  }
305  bfcur++;
306  }
307  }
308 
309  // now y component basis functions
310  for (int j=0;j<openBasis_.getCardinality();j++)
311  {
312  for (int i=0;i<closedBasis_.getCardinality();i++)
313  {
314  for (int l=0;l<dim0;l++)
315  {
316  outputValues(bfcur,l) = closedBasisDerivsXPts(i,l,0) * openBasisValsYPts(j,l);
317  }
318  bfcur++;
319  }
320  }
321  }
322  break;
323  case OPERATOR_DIV:
324  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
325  ">>> ERROR (Basis_HCURL_QUAD_In_FEM): DIV is invalid operator for HCURL Basis Functions");
326  break;
327 
328  case OPERATOR_GRAD:
329  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_GRAD), std::invalid_argument,
330  ">>> ERROR (Basis_HCURL_QUAD_In_FEM): GRAD is invalid operator for HCURL Basis Functions");
331  break;
332 
333  case OPERATOR_D1:
334  case OPERATOR_D2:
335  case OPERATOR_D3:
336  case OPERATOR_D4:
337  case OPERATOR_D5:
338  case OPERATOR_D6:
339  case OPERATOR_D7:
340  case OPERATOR_D8:
341  case OPERATOR_D9:
342  case OPERATOR_D10:
343  TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType == OPERATOR_D1) ||
344  (operatorType == OPERATOR_D2) ||
345  (operatorType == OPERATOR_D3) ||
346  (operatorType == OPERATOR_D4) ||
347  (operatorType == OPERATOR_D5) ||
348  (operatorType == OPERATOR_D6) ||
349  (operatorType == OPERATOR_D7) ||
350  (operatorType == OPERATOR_D8) ||
351  (operatorType == OPERATOR_D9) ||
352  (operatorType == OPERATOR_D10) ),
353  std::invalid_argument,
354  ">>> ERROR (Basis_HCURL_QUAD_In_FEM): Invalid operator type");
355  break;
356 
357  default:
358  TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType != OPERATOR_VALUE) &&
359  (operatorType != OPERATOR_GRAD) &&
360  (operatorType != OPERATOR_CURL) &&
361  (operatorType != OPERATOR_CURL) &&
362  (operatorType != OPERATOR_D1) &&
363  (operatorType != OPERATOR_D2) &&
364  (operatorType != OPERATOR_D3) &&
365  (operatorType != OPERATOR_D4) &&
366  (operatorType != OPERATOR_D5) &&
367  (operatorType != OPERATOR_D6) &&
368  (operatorType != OPERATOR_D7) &&
369  (operatorType != OPERATOR_D8) &&
370  (operatorType != OPERATOR_D9) &&
371  (operatorType != OPERATOR_D10) ),
372  std::invalid_argument,
373  ">>> ERROR (Basis_HCURL_QUAD_In_FEM): Invalid operator type");
374  }
375  }
376 
377 
378 
379  template<class Scalar, class ArrayScalar>
381  const ArrayScalar & inputPoints,
382  const ArrayScalar & cellVertices,
383  const EOperator operatorType) const {
384  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
385  ">>> ERROR (Basis_HCURL_QUAD_In_FEM): FEM Basis calling an FVD member function");
386  }
387 
388  template<class Scalar, class ArrayScalar>
390  {
391  // x-component basis functions
392  int cur = 0;
393 
394  for (int j=0;j<closedPts_.dimension(0);j++)
395  {
396  for (int i=0;i<openPts_.dimension(0);i++)
397  {
398  DofCoords(cur,0) = openPts_(i,0);
399  DofCoords(cur,1) = closedPts_(j,0);
400  cur++;
401  }
402  }
403 
404  // y-component basis functions
405  for (int j=0;j<openPts_.dimension(0);j++)
406  {
407  for (int i=0;i<closedPts_.dimension(0);i++)
408  {
409  DofCoords(cur,0) = closedPts_(i,0);
410  DofCoords(cur,1) = openPts_(j,0);
411  cur++;
412  }
413  }
414 
415  return;
416  }
417 
418 
419 
420 }// namespace Intrepid
421 
422 #if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
423 #ifdef __GNUC__
424 #warning "The Intrepid package is deprecated"
425 #endif
426 #endif
427 
EBasis basisType_
Type of the basis.
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
virtual void getDofCoords(ArrayScalar &DofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on a reference cell; defined for interp...
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Quadrilateral cell.
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package http://trilinos...
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
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. ...
Basis_HCURL_QUAD_In_FEM(int order, const ArrayScalar &ptsClosed, const ArrayScalar &ptsOpen)
Constructor.