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
49namespace 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>
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
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.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Quadrilateral cell.
Basis_HDIV_QUAD_In_FEM(int order, const ArrayScalar &ptsClosed, const ArrayScalar &ptsOpen)
Constructor.
virtual void getDofCoords(ArrayScalar &DofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on a reference cell; defined for interp...
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
bool basisTagsAreSet_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
int basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom.
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
EBasis basisType_
Type of the basis.
int basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package http://trilinos....
static void lineProduct2d(const int dim0, const int entity0, const int dim1, const int entity1, int &resultdim, int &resultentity)