|
Intrepid
|
00001 #include "Teuchos_Array.hpp" 00002 #include "Teuchos_RCP.hpp" 00003 #include "Teuchos_oblackholestream.hpp" 00004 #include "Teuchos_GlobalMPISession.hpp" 00005 #include "Intrepid_TensorProductSpaceTools.hpp" 00006 #include "Intrepid_HGRAD_QUAD_Cn_FEM.hpp" 00007 #include "Intrepid_HGRAD_HEX_Cn_FEM.hpp" 00008 #include "Intrepid_CubaturePolylib.hpp" 00009 #include "Intrepid_Utils.hpp" 00010 #include "Intrepid_Types.hpp" 00011 00012 00013 using Teuchos::Array; 00014 using Intrepid::FieldContainer; 00015 using Intrepid::Basis; 00016 using Intrepid::TensorBasis; 00017 00018 #define INTREPID_TEST_COMMAND( S ) \ 00019 { \ 00020 try { \ 00021 S ; \ 00022 } \ 00023 catch (std::logic_error err) { \ 00024 *outStream << "Expected Error ----------------------------------------------------------------\n"; \ 00025 *outStream << err.what() << '\n'; \ 00026 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \ 00027 }; \ 00028 } 00029 00030 00031 int main( int argc , char **argv ) 00032 { 00033 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00034 00035 // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided. 00036 int iprint = argc - 1; 00037 00038 Teuchos::RCP<std::ostream> outStream; 00039 Teuchos::oblackholestream bhs; // outputs nothing 00040 00041 if (iprint > 0) 00042 outStream = Teuchos::rcp(&std::cout, false); 00043 else 00044 outStream = Teuchos::rcp(&bhs, false); 00045 00046 // Save the format state of the original std::cout. 00047 Teuchos::oblackholestream oldFormatState; 00048 oldFormatState.copyfmt(std::cout); 00049 00050 00051 *outStream \ 00052 << "===============================================================================\n" \ 00053 << "| |\n" \ 00054 << "| Unit Test TensorProductSpace Tools |\n" \ 00055 << "| |\n" \ 00056 << "| Tests sum-factored polynomial evaluation and integration |\n" \ 00057 << "| |\n" \ 00058 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \ 00059 << "| Denis Ridzal (dridzal@sandia.gov) or |\n" \ 00060 << "| Robert Kirby (robert.c.kirby@ttu.edu) |\n" \ 00061 << "| |\n" \ 00062 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00063 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00064 << "| |\n" \ 00065 << "===============================================================================\n"; 00066 00067 00068 00069 int errorFlag = 0; 00070 00071 Array<RCP<TensorBasis<double,FieldContainer<double> > > > basesByDim(4); 00072 Array<RCP<FieldContainer<double> > > cubPtsByDim(4); 00073 00074 Intrepid::CubaturePolylib<double> cpl(2,Intrepid::PL_GAUSS_LOBATTO); 00075 00076 FieldContainer<double> cubPoints( cpl.getNumPoints() ,1 ); 00077 FieldContainer<double> cubWeights( cpl.getNumPoints() ); 00078 00079 cpl.getCubature( cubPoints, cubWeights ); 00080 00081 basesByDim[2] = Teuchos::rcp( new Intrepid::Basis_HGRAD_QUAD_Cn_FEM<double,FieldContainer<double> >( 2 , Intrepid::POINTTYPE_SPECTRAL ) ); 00082 basesByDim[3] = Teuchos::rcp( new Intrepid::Basis_HGRAD_HEX_Cn_FEM<double,FieldContainer<double> >( 2 , Intrepid::POINTTYPE_SPECTRAL ) ); 00083 00084 00085 // get points 00086 FieldContainer<double> quadPts( cpl.getNumPoints() * cpl.getNumPoints() , 2 ); 00087 for (int j=0;j<cpl.getNumPoints();j++) 00088 { 00089 for (int i=0;i<cpl.getNumPoints();i++) 00090 { 00091 int index = j*cpl.getNumPoints() + i; 00092 quadPts(index,0) = cubPoints(i,0); 00093 quadPts(index,1) = cubPoints(j,0); 00094 } 00095 } 00096 00097 FieldContainer<double> cubPts( cpl.getNumPoints() * cpl.getNumPoints() * cpl.getNumPoints() , 3 ); 00098 for (int k=0;k<cpl.getNumPoints();k++) 00099 { 00100 for (int j=0;j<cpl.getNumPoints();j++) 00101 { 00102 for (int i=0;i<cpl.getNumPoints();i++) 00103 { 00104 int index = k* cpl.getNumPoints() * cpl.getNumPoints() + j*cpl.getNumPoints() + i; 00105 cubPts(index,0) = cubPoints(i,0); 00106 cubPts(index,1) = cubPoints(j,0); 00107 cubPts(index,2) = cubPoints(k,0); 00108 } 00109 } 00110 } 00111 00112 cubPtsByDim[2] = Teuchos::rcp( &quadPts , false ); 00113 cubPtsByDim[3] = Teuchos::rcp( &cubPts , false ); 00114 00115 int space_dim = 2; 00116 00117 Array<Array<RCP<Basis<double,FieldContainer<double> > > > > &bases = basesByDim[space_dim]->getBases(); 00118 00119 FieldContainer<double> coeff(1,1,basesByDim[space_dim]->getCardinality()); 00120 00121 00122 00123 Array<RCP<FieldContainer<double> > > pts( space_dim ); 00124 pts[0] = Teuchos::rcp( &cubPoints, false ); 00125 for (int i=1;i<space_dim;i++) 00126 { 00127 pts[i] = pts[0]; 00128 } 00129 00130 Array<RCP<FieldContainer<double> > > wts(space_dim); 00131 wts[0] = Teuchos::rcp( &cubWeights , false ); 00132 for (int i=1;i<space_dim;i++) 00133 { 00134 wts[i] = wts[0]; 00135 } 00136 00137 FieldContainer<double> Phix(bases[0][0]->getCardinality(), 00138 cpl.getNumPoints() ); 00139 FieldContainer<double> Phiy(bases[0][1]->getCardinality(), 00140 cpl.getNumPoints() ); 00141 FieldContainer<double> DPhix(bases[0][0]->getCardinality(), 00142 cpl.getNumPoints(), 1 ); 00143 FieldContainer<double> DPhiy(bases[0][1]->getCardinality(), 00144 cpl.getNumPoints(), 1 ); 00145 00146 bases[0][0]->getValues( Phix , cubPoints, Intrepid::OPERATOR_VALUE ); 00147 bases[0][1]->getValues( Phiy , cubPoints, Intrepid::OPERATOR_VALUE ); 00148 bases[0][0]->getValues( DPhix , cubPoints, Intrepid::OPERATOR_D1 ); 00149 bases[0][1]->getValues( DPhiy , cubPoints, Intrepid::OPERATOR_D1 ); 00150 00151 Array<RCP<FieldContainer<double> > > basisVals(2); 00152 basisVals[0] = Teuchos::rcp( &Phix , false ); 00153 basisVals[1] = Teuchos::rcp( &Phiy , false ); 00154 00155 Array<RCP<FieldContainer<double> > > basisDVals(2); 00156 basisDVals[0] = Teuchos::rcp( &DPhix , false ); 00157 basisDVals[1] = Teuchos::rcp( &DPhiy , false ); 00158 00159 FieldContainer<double> vals(1,1,pts[0]->size() * pts[1]->size() ); 00160 00161 // first basis function is the polynomial. 00162 coeff(0,0,0) = 1.0; 00163 00164 Intrepid::TensorProductSpaceTools::evaluate<double,FieldContainer<double>,FieldContainer<double>,FieldContainer<double> >( vals , coeff , basisVals ); 00165 00166 FieldContainer<double> grads( 1 , 1, pts[0]->size() * pts[1]->size() , 2 ); 00167 00168 Intrepid::TensorProductSpaceTools::evaluateGradient<double,FieldContainer<double>,FieldContainer<double>,FieldContainer<double> >( grads , 00169 coeff , 00170 basisVals , 00171 basisDVals ); 00172 00173 // confirm by comparing to actual gradients 00174 FieldContainer<double> fullVals(basesByDim[space_dim]->getCardinality(), 00175 basesByDim[space_dim]->getCardinality()); 00176 FieldContainer<double> fullGrads(basesByDim[space_dim]->getCardinality(), 00177 basesByDim[space_dim]->getCardinality(), 00178 space_dim ); 00179 00180 00181 basesByDim[space_dim]->getValues( fullVals , 00182 quadPts , 00183 Intrepid::OPERATOR_VALUE ); 00184 basesByDim[space_dim]->getValues( fullGrads , 00185 quadPts , 00186 Intrepid::OPERATOR_GRAD ); 00187 00188 for (int i=0;i<fullVals.dimension(1);i++) 00189 { 00190 if (std::abs( fullVals(0,i) - vals(0,0,i) ) > Intrepid::INTREPID_TOL ) 00191 { 00192 errorFlag++; 00193 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00194 00195 // Output the multi-index of the value where the error is: 00196 *outStream << " Evaluating first bf At multi-index { "; 00197 *outStream << i; 00198 *outStream << "} brute force value: " << fullVals(0,i) 00199 << " but tensor-product value: " << vals(0,i) << "\n"; 00200 *outStream << "Difference: " << std::abs( fullVals(0,i) - vals(0,i) ) << "\n"; 00201 } 00202 } 00203 00204 for (int i=0;i<fullGrads.dimension(1);i++) 00205 { 00206 for (int j=0;j<fullGrads.dimension(2);j++) 00207 { 00208 if (std::abs( fullGrads(0,i,j) - grads(0,0,i,j) ) > Intrepid::INTREPID_TOL ) 00209 { 00210 errorFlag++; 00211 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00212 00213 // Output the multi-index of the value where the error is: 00214 *outStream << " Evaluating first bf At multi-index { "; 00215 *outStream << i << " " << j; 00216 *outStream << "} brute force value: " << fullGrads(0,i,j) 00217 << " but tensor-product value: " << grads(0,0,i,j) << "\n"; 00218 *outStream << "Difference: " << std::abs( fullGrads(0,i,j) - grads(0,0,i,j) ) << "\n"; 00219 } 00220 } 00221 } 00222 00223 00224 // now test moments. 00225 // I've already evaluated the first basis function at the quadrature points. 00226 // why not use it? 00227 00228 FieldContainer<double> momentsNaive(1,basesByDim[2]->getCardinality()); 00229 for (int i=0;i<basesByDim[2]->getCardinality();i++) 00230 { 00231 momentsNaive(0,i) = 0.0; 00232 for (int qpty=0;qpty<cubPoints.dimension(0);qpty++) 00233 { 00234 for (int qptx=0;qptx<cubPoints.dimension(0);qptx++) 00235 { 00236 momentsNaive(0,i) += cubWeights(qpty) * cubWeights(qptx) * 00237 vals( 0, 0, qpty*cubPoints.dimension(0)+qptx ) 00238 * fullVals(i,qpty*cubPoints.dimension(0)+qptx); 00239 } 00240 } 00241 } 00242 00243 FieldContainer<double> momentsClever(1,1,basesByDim[space_dim]->getCardinality()); 00244 Intrepid::TensorProductSpaceTools::moments<double,FieldContainer<double>,FieldContainer<double>,FieldContainer<double>,FieldContainer<double> >( momentsClever , 00245 vals , 00246 basisVals , 00247 wts ); 00248 for (int j=0;j<momentsClever.dimension(0);j++) 00249 { 00250 if (std::abs( momentsClever(0,0,j) - momentsNaive(0,j) ) > Intrepid::INTREPID_TOL ) 00251 { 00252 errorFlag++; 00253 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00254 00255 // Output the multi-index of the value where the error is: 00256 *outStream << " At multi-index { "; 00257 *outStream << " " << j; 00258 *outStream << "} brute force value: " << momentsNaive(0,j) 00259 << " but sum-factored value: " << momentsClever(0,0,j) << "\n"; 00260 *outStream << "Difference: " << std::abs( momentsNaive(0,j) - momentsClever(0,0,j) ) << "\n"; 00261 } 00262 } 00263 00264 if (errorFlag != 0) 00265 { 00266 std::cout << "End Result: TEST FAILED\n"; 00267 } 00268 else 00269 { 00270 std::cout << "End Result: TEST PASSED\n"; 00271 } 00272 00273 // reset format state of std::cout 00274 std::cout.copyfmt(oldFormatState); 00275 00276 return errorFlag; 00277 00278 }
1.7.6.1