本文整理汇总了C++中VarFactoryPtr类的典型用法代码示例。如果您正苦于以下问题:C++ VarFactoryPtr类的具体用法?C++ VarFactoryPtr怎么用?C++ VarFactoryPtr使用的例子?那么恭喜您, 这里精选的类代码示例或许可以为您提供帮助。
在下文中一共展示了VarFactoryPtr类的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: quadPoints
void RHSTests::setup()
{
FieldContainer<double> quadPoints(4,2);
quadPoints(0,0) = 0.0; // x1
quadPoints(0,1) = 0.0; // y1
quadPoints(1,0) = 1.0;
quadPoints(1,1) = 0.0;
quadPoints(2,0) = 1.0;
quadPoints(2,1) = 1.0;
quadPoints(3,0) = 0.0;
quadPoints(3,1) = 1.0;
int H1Order = 3;
int delta_p = 3; // for test functions
int horizontalCells = 2;
int verticalCells = 2;
double eps = 1.0; // not really testing for sharp gradients right now--just want to see if things basically work
double beta_x = 1.0;
double beta_y = 1.0;
BFPtr confusionBF = ConfusionBilinearForm::confusionBF(eps,beta_x,beta_y);
Teuchos::RCP<ConfusionProblemLegacy> confusionProblem = Teuchos::rcp( new ConfusionProblemLegacy(confusionBF, beta_x, beta_y) );
_rhs = confusionProblem;
_mesh = MeshFactory::buildQuadMesh(quadPoints, horizontalCells, verticalCells, confusionBF, H1Order, H1Order+delta_p);
_mesh->setUsePatchBasis(false);
VarFactoryPtr varFactory = confusionBF->varFactory(); // Create test IDs that match the enum in ConfusionBilinearForm
_tau = varFactory->testVar(ConfusionBilinearForm::S_TAU,HDIV);
_v = varFactory->testVar(ConfusionBilinearForm::S_V,HGRAD);
_rhsEasy = RHS::rhs();
_rhsEasy->addTerm( _v );
}
开发者ID:vijaysm,项目名称:Camellia,代码行数:35,代码来源:RHSTests.cpp
示例2: phi
void PoissonExactSolution::setUseSinglePointBCForPHI(bool useSinglePointBCForPhi, IndexType vertexIndexForZeroValue)
{
FunctionPtr phi_exact = phi();
VarFactoryPtr vf = _bf->varFactory();
VarPtr psi_hat_n = vf->fluxVar(PoissonBilinearForm::S_PSI_HAT_N);
VarPtr q = vf->testVar(PoissonBilinearForm::S_Q, HGRAD);
VarPtr phi = vf->fieldVar(PoissonBilinearForm::S_PHI);
SpatialFilterPtr wholeBoundary = SpatialFilter::allSpace();
FunctionPtr n = Function::normal();
FunctionPtr psi_n_exact = phi_exact->grad() * n;
_bc = BC::bc();
_bc->addDirichlet(psi_hat_n, wholeBoundary, psi_n_exact);
if (!useSinglePointBCForPhi)
{
_bc->addZeroMeanConstraint(phi);
}
else
{
std::vector<double> point = getPointForBCImposition();
double value = Function::evaluate(phi_exact, point[0], point[1]);
// cout << "PoissonExactSolution: imposing phi = " << value << " at (" << point[0] << ", " << point[1] << ")\n";
_bc->addSpatialPointBC(phi->ID(), value, point);
}
}
开发者ID:vijaysm,项目名称:Camellia,代码行数:30,代码来源:PoissonExactSolution.cpp
示例3: tau
VarPtr PoissonFormulation::tau()
{
VarFactoryPtr vf = _poissonBF->varFactory();
if (_spaceDim > 1)
return vf->testVar(S_TAU, HDIV);
else
return vf->testVar(S_TAU, HGRAD);
}
开发者ID:CamelliaDPG,项目名称:Camellia,代码行数:8,代码来源:PoissonFormulation.cpp
示例4: bf
BFPtr TestBilinearFormDx::bf()
{
VarFactoryPtr vf = VarFactory::varFactory();
VarPtr u = vf->fieldVar("u",HGRAD);
VarPtr v = vf->testVar("v",HGRAD);
BFPtr bf = BF::bf(vf);
bf->addTerm(u->dx(), v->dx());
return bf;
}
开发者ID:CamelliaDPG,项目名称:Camellia,代码行数:10,代码来源:TestBilinearFormDx.cpp
示例5: testLinearTermEvaluation
bool LinearTermTests::testLinearTermEvaluation()
{
bool success = true;
double eps = .1;
FunctionPtr one = Function::constant(1.0);
vector<double> e1,e2;
e1.push_back(1.0);
e1.push_back(0.0);
e2.push_back(0.0);
e2.push_back(1.0);
// define test variables
VarFactoryPtr varFactory = VarFactory::varFactory();
VarPtr tau = varFactory->testVar("\\tau", HDIV);
VarPtr v = varFactory->testVar("v", HGRAD);
// define a couple LinearTerms
LinearTermPtr vVecLT = Teuchos::rcp(new LinearTerm);
LinearTermPtr tauVecLT = Teuchos::rcp(new LinearTerm);
vVecLT->addTerm(sqrt(eps)*v->grad());
tauVecLT->addTerm((1/sqrt(eps))*tau);
//////////////////// evaluate LinearTerms /////////////////
map<int,FunctionPtr> errRepMap;
errRepMap[v->ID()] = one;
errRepMap[tau->ID()] = one*e1+one*e2; // vector valued fxn (1,1)
FunctionPtr errTau = tauVecLT->evaluate(errRepMap,false);
FunctionPtr errV = vVecLT->evaluate(errRepMap,false);
try
{
bool xTauZero = errTau->x()->isZero();
bool yTauZero = errTau->y()->isZero();
bool xVZero = errV->dx()->isZero();
bool yVZero = errV->dy()->isZero();
}
catch (...)
{
cout << "testLinearTermEvaluation: Caught exception.\n";
success = false;
}
/*
FunctionPtr xErr = (errTau->x())*(errTau->x()) + (errV->dx())*(errV->dx());
FunctionPtr yErr = (errTau->y())*(errTau->y()) + (errV->dy())*(errV->dy());
double xErrVal = xErr->integrate(mesh,15,true);
*/
// if we don't crash, return success
return success;
}
开发者ID:CamelliaDPG,项目名称:Camellia,代码行数:53,代码来源:LinearTermTests.cpp
示例6: testLobattoValues
bool LobattoBasisTests::testLobattoValues()
{
bool success = true;
FunctionPtr x = Function::xn(1);
vector< FunctionPtr > lobattoFunctionsExpected;
// Pavel Solin's first two Lobatto functions:
// lobattoFunctionsExpected.push_back( (1 - x) / 2 );
// lobattoFunctionsExpected.push_back( (1 + x) / 2 );
// Demkowicz's:
lobattoFunctionsExpected.push_back( Function::constant(1.0) );
lobattoFunctionsExpected.push_back( x );
lobattoFunctionsExpected.push_back( (x*x - 1) / 2);
lobattoFunctionsExpected.push_back( (x*x - 1) * x / 2);
lobattoFunctionsExpected.push_back( (x*x - 1) * (5 * x * x - 1) / 8);
lobattoFunctionsExpected.push_back( (x*x - 1) * (7 * x * x - 3) * x / 8);
vector< FunctionPtr > lobattoFunctions;
bool conformingFalse = false;
int n_max = 4;
for (int n=0; n<=n_max; n++)
{
lobattoFunctions.push_back( Teuchos::rcp( new LobattoFunction<>(n,conformingFalse,false) ) );
}
VarFactoryPtr varFactory = VarFactory::varFactory();
varFactory->testVar("v", HGRAD);
varFactory->fieldVar("u", L2);
BFPtr bf = Teuchos::rcp( new BF(varFactory) );
MeshPtr mesh = MeshFactory::quadMesh(bf, n_max);
BasisCachePtr basisCache = BasisCache::basisCacheForCell(mesh, 0);
double tol = 1e-12;
for (int n=0; n<=n_max; n++)
{
if (! lobattoFunctions[n]->equals(lobattoFunctionsExpected[n], basisCache, tol) )
{
cout << "Lobatto function " << n << " does not match expected.\n";
success = false;
}
}
return success;
}
开发者ID:CamelliaDPG,项目名称:Camellia,代码行数:49,代码来源:LobattoBasisTests.cpp
示例7: testLobattoDerivativeValues
bool LobattoBasisTests::testLobattoDerivativeValues()
{
bool success = true;
FunctionPtr x = Function::xn(1);
vector< FunctionPtr > legendreFunctions; // manually specified
vector< FunctionPtr > lobattoDerivatives;
FunctionPtr one = Function::constant(1);
legendreFunctions.push_back(one);
legendreFunctions.push_back(x);
legendreFunctions.push_back(0.5 * (3 * x * x - 1) );
legendreFunctions.push_back( (5 * x * x * x - 3 * x) / 2);
legendreFunctions.push_back( (1.0 / 8.0) * (35 * x * x * x * x - 30 * x * x + 3) );
legendreFunctions.push_back( (1.0 / 8.0) * (63 * x * x * x * x * x - 70 * x * x * x + 15 * x) );
bool conformingFalse = false;
int n_max = 4;
for (int n=0; n<=n_max+1; n++)
{
lobattoDerivatives.push_back( Teuchos::rcp( new LobattoFunction<>(n,conformingFalse,true) ) );
}
VarFactoryPtr varFactory = VarFactory::varFactory();
varFactory->testVar("v", HGRAD);
varFactory->fieldVar("u", L2);
BFPtr bf = Teuchos::rcp( new BF(varFactory) );
MeshPtr mesh = MeshFactory::quadMesh(bf, n_max);
BasisCachePtr basisCache = BasisCache::basisCacheForCell(mesh, 0);
double tol = 1e-8;
for (int n=0; n<=n_max; n++)
{
if (! legendreFunctions[n]->equals(lobattoDerivatives[n+1], basisCache, tol) )
{
cout << "Legendre function " << n << " != Lobatto function " << n + 1 << " derivative.\n";
cout << "L_" << n << "(0.5) = " << Function::evaluate(legendreFunctions[n],0.5) << endl;
cout << "l'_" << n+1 << "(0.5) = " << Function::evaluate(lobattoDerivatives[n+1],0.5) << endl;
success = false;
}
}
return success;
}
开发者ID:CamelliaDPG,项目名称:Camellia,代码行数:46,代码来源:LobattoBasisTests.cpp
示例8: TEUCHOS_TEST_FOR_EXCEPTION
VarPtr PressurelessStokesFormulation::v(int i)
{
if (i > _spaceDim)
{
TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "i must be less than or equal to _spaceDim");
}
VarFactoryPtr vf = _stokesBF->varFactory();
switch (i)
{
case 1:
return vf->testVar(S_V1, HGRAD);
case 2:
return vf->testVar(S_V2, HGRAD);
case 3:
return vf->testVar(S_V3, HGRAD);
}
TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "unhandled i value");
}
开发者ID:CamelliaDPG,项目名称:Camellia,代码行数:18,代码来源:PressurelessStokesFormulation.cpp
示例9: LegendreFunction
bool LobattoBasisTests::testLegendreValues()
{
bool success = true;
FunctionPtr x = Function::xn(1);
vector< FunctionPtr > legendreFunctionsExpected;
legendreFunctionsExpected.push_back( Function::constant(1.0) );
legendreFunctionsExpected.push_back( x );
legendreFunctionsExpected.push_back( (3 * x*x - 1) / 2);
legendreFunctionsExpected.push_back( (5 * x*x*x - 3*x) / 2);
legendreFunctionsExpected.push_back( (35 * x * x * x * x - 30 * x * x + 3) / 8);
legendreFunctionsExpected.push_back( (63 * x * x * x * x * x - 70 * x * x * x + 15 * x) / 8);
vector< FunctionPtr > legendreFunctions;
int n_max = 4;
for (int n=0; n<=n_max; n++)
{
legendreFunctions.push_back( Teuchos::rcp( new LegendreFunction(n) ) );
}
VarFactoryPtr varFactory = VarFactory::varFactory();
varFactory->testVar("v", HGRAD);
varFactory->fieldVar("u", L2);
BFPtr bf = Teuchos::rcp( new BF(varFactory) );
MeshPtr mesh = MeshFactory::quadMesh(bf, n_max);
BasisCachePtr basisCache = BasisCache::basisCacheForCell(mesh, 0);
double tol = 1e-8; // relax to make sure that failure isn't just roundoff
for (int n=0; n<=n_max; n++)
{
if (! legendreFunctions[n]->equals(legendreFunctionsExpected[n], basisCache, tol) )
{
cout << "Legendre function " << n << " does not match expected.\n";
success = false;
}
}
return success;
}
开发者ID:CamelliaDPG,项目名称:Camellia,代码行数:42,代码来源:LobattoBasisTests.cpp
示例10:
ConfusionManufacturedSolution::ConfusionManufacturedSolution(double epsilon, double beta_x, double beta_y)
{
_epsilon = epsilon;
_beta_x = beta_x;
_beta_y = beta_y;
// set the class variables from ExactSolution:
// _bc = Teuchos::rcp(this,false); // false: don't let the RCP own the memory
// _rhs = Teuchos::rcp(this,false);
BFPtr bf = ConfusionBilinearForm::confusionBF(epsilon,beta_x,beta_y);
_bilinearForm = bf;
VarFactoryPtr vf = bf->varFactory();
VarPtr u = vf->fieldVar(ConfusionBilinearForm::S_U);
VarPtr sigma1 = vf->fieldVar(ConfusionBilinearForm::S_SIGMA_1);
VarPtr sigma2 = vf->fieldVar(ConfusionBilinearForm::S_SIGMA_2);
_u_hat = vf->traceVar(ConfusionBilinearForm::S_U_HAT);
_beta_n_u_minus_sigma_hat = vf->fluxVar(ConfusionBilinearForm::S_BETA_N_U_MINUS_SIGMA_HAT);
_v = vf->testVar(ConfusionBilinearForm::S_V, HGRAD);
FunctionPtr u_exact = this->u();
FunctionPtr sigma_exact = epsilon * u_exact->grad();
FunctionPtr u_exact_laplacian = u_exact->dx()->dx() + u_exact->dy()->dy();
_rhs = RHS::rhs();
FunctionPtr f = - _epsilon * u_exact_laplacian + _beta_x * u_exact->dx() + _beta_y * u_exact->dy();
_rhs->addTerm( f * _v );
_bc = BC::bc();
_bc->addDirichlet(_u_hat, SpatialFilter::allSpace(), u_exact);
FunctionPtr beta = Function::vectorize(Function::constant(_beta_x), Function::constant(_beta_y));
FunctionPtr n = Function::normal();
FunctionPtr one_skeleton = Function::meshSkeletonCharacteristic(); // allows restriction to skeleton
FunctionPtr sigma_flux_exact = beta * ( n * u_exact - sigma_exact * one_skeleton);
this->setSolutionFunction(u, u_exact);
this->setSolutionFunction(sigma1, sigma_exact->x());
this->setSolutionFunction(sigma2, sigma_exact->y());
this->setSolutionFunction(_u_hat, u_exact);
this->setSolutionFunction(_beta_n_u_minus_sigma_hat, sigma_flux_exact);
}
开发者ID:CamelliaDPG,项目名称:Camellia,代码行数:48,代码来源:ConfusionManufacturedSolution.cpp
示例11: setup
void LinearTermTests::setup()
{
// VarPtr v1, v2, v3; // HGRAD members (test variables)
// VarPtr q1, q2, q3; // HDIV members (test variables)
// VarPtr u1, u2, u3; // L2 members (trial variables)
// VarPtr u1_hat, u2_hat; // trace variables
// VarPtr u3_hat_n; // flux variable
//
// FunctionPtr sine_x;
sine_x = Teuchos::rcp( new Sine_x );
cos_y = Teuchos::rcp( new Cosine_y );
VarFactoryPtr varFactory = VarFactory::varFactory();
q1 = varFactory->testVar("q_1", HDIV);
q2 = varFactory->testVar("q_2", HDIV);
q3 = varFactory->testVar("q_3", HDIV);
v1 = varFactory->testVar("v_1", HGRAD);
v2 = varFactory->testVar("v_2", HGRAD);
v3 = varFactory->testVar("v_3", HGRAD);
u1 = varFactory->fieldVar("u_1", HGRAD);
u2 = varFactory->fieldVar("u_2", HGRAD);
u3 = varFactory->fieldVar("u_3", HGRAD);
u1_hat = varFactory->traceVar("\\widehat{u}_1");
u2_hat = varFactory->traceVar("\\widehat{u}_2");
u3_hat_n = varFactory->fluxVar("\\widehat{u}_3n");
bf = Teuchos::rcp(new BF(varFactory)); // made-up bf for Mesh + previous solution tests
bf->addTerm(u1_hat, q1->dot_normal());
bf->addTerm(u1, q1->x());
bf->addTerm(u2, q1->y());
bf->addTerm(u3_hat_n, v1);
bf->addTerm(u3, v1);
// DofOrderingFactory discreteSpaceFactory(bf);
int polyOrder = 3, testToAdd = 2;
Teuchos::RCP<shards::CellTopology> quadTopoPtr;
// quadTopoPtr = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() ));
// define nodes for mesh
FieldContainer<double> quadPoints(4,2);
quadPoints(0,0) = -1.0; // x1
quadPoints(0,1) = -1.0; // y1
quadPoints(1,0) = 1.0;
quadPoints(1,1) = -1.0;
quadPoints(2,0) = 1.0;
quadPoints(2,1) = 1.0;
quadPoints(3,0) = -1.0;
quadPoints(3,1) = 1.0;
int horizontalElements = 2, verticalElements = 2;
mesh = MeshFactory::buildQuadMesh(quadPoints, horizontalElements, verticalElements, bf, polyOrder+1, polyOrder+1+testToAdd);
ElementTypePtr elemType = mesh->getElement(0)->elementType();
trialOrder = elemType->trialOrderPtr;
testOrder = elemType->testOrderPtr;
basisCache = Teuchos::rcp(new BasisCache(elemType, mesh));
vector<GlobalIndexType> cellIDs;
cellIDs.push_back(0);
cellIDs.push_back(1);
cellIDs.push_back(2);
cellIDs.push_back(3);
bool createSideCacheToo = true;
basisCache->setPhysicalCellNodes(mesh->physicalCellNodesGlobal(elemType), cellIDs, createSideCacheToo);
}
开发者ID:CamelliaDPG,项目名称:Camellia,代码行数:77,代码来源:LinearTermTests.cpp
示例12: testIntegrateMixedBasis
bool LinearTermTests::testIntegrateMixedBasis()
{
bool success = true;
//////////////////// DECLARE VARIABLES ///////////////////////
// define test variables
VarFactoryPtr varFactory = VarFactory::varFactory();
VarPtr v = varFactory->testVar("v", HGRAD);
// define trial variables
VarPtr beta_n_u_hat = varFactory->fluxVar("\\widehat{\\beta \\cdot n }");
VarPtr u = varFactory->fieldVar("u");
vector<double> beta;
beta.push_back(1.0);
beta.push_back(1.0);
//////////////////// DEFINE BILINEAR FORM/Mesh ///////////////////////
BFPtr convectionBF = Teuchos::rcp( new BF(varFactory) );
// v terms:
convectionBF->addTerm( -u, beta * v->grad() );
convectionBF->addTerm( beta_n_u_hat, v);
convectionBF->addTerm( u, v);
// build CONSTANT SINGLE ELEMENT MESH
int order = 0;
int H1Order = order+1;
int pToAdd = 1;
int nCells = 2; // along a side
// create a pointer to a new mesh:
Teuchos::RCP<Mesh> mesh = MeshUtilities::buildUnitQuadMesh(nCells,convectionBF, H1Order, H1Order+pToAdd);
ElementTypePtr elemType = mesh->getElement(0)->elementType();
BasisCachePtr basisCache = Teuchos::rcp(new BasisCache(elemType, mesh));
vector<GlobalIndexType> cellIDs;
vector< ElementPtr > allElems = mesh->activeElements();
vector< ElementPtr >::iterator elemIt;
for (elemIt=allElems.begin(); elemIt!=allElems.end(); elemIt++)
{
cellIDs.push_back((*elemIt)->cellID());
}
bool createSideCacheToo = true;
basisCache->setPhysicalCellNodes(mesh->physicalCellNodesGlobal(elemType), cellIDs, createSideCacheToo);
int numTrialDofs = elemType->trialOrderPtr->totalDofs();
int numCells = mesh->numActiveElements();
double areaPerCell = 1.0 / numCells;
FieldContainer<double> integrals(numCells,numTrialDofs);
FieldContainer<double> expectedIntegrals(numCells,numTrialDofs);
double sidelengthOfCell = 1.0 / nCells;
DofOrderingPtr trialOrdering = elemType->trialOrderPtr;
int dofForField = trialOrdering->getDofIndex(u->ID(), 0);
vector<int> dofsForFlux;
const vector<int>* sidesForFlux = &trialOrdering->getSidesForVarID(beta_n_u_hat->ID());
for (vector<int>::const_iterator sideIt = sidesForFlux->begin(); sideIt != sidesForFlux->end(); sideIt++)
{
int sideIndex = *sideIt;
dofsForFlux.push_back(trialOrdering->getDofIndex(beta_n_u_hat->ID(), 0, sideIndex));
}
for (int cellIndex = 0; cellIndex < numCells; cellIndex++)
{
expectedIntegrals(cellIndex, dofForField) = areaPerCell;
for (vector<int>::iterator dofIt = dofsForFlux.begin(); dofIt != dofsForFlux.end(); dofIt++)
{
int fluxDofIndex = *dofIt;
expectedIntegrals(cellIndex, fluxDofIndex) = sidelengthOfCell;
}
}
// cout << "expectedIntegrals:\n" << expectedIntegrals;
// setup: with constant degrees of freedom, we expect that the integral of int_dK (flux) + int_K (field) will be ones for each degree of freedom, assuming the basis functions for these constants field/flux variables are just C = 1.0.
//
//On a unit square, int_K (constant) = 1.0, and int_dK (u_i) = 1, for i = 0,...,3.
LinearTermPtr lt = 1.0 * beta_n_u_hat;
LinearTermPtr field = 1.0 * u;
lt->addTerm(field,true);
lt->integrate(integrals, elemType->trialOrderPtr, basisCache);
double tol = 1e-12;
double maxDiff;
success = TestSuite::fcsAgree(integrals,expectedIntegrals,tol,maxDiff);
if (success==false)
{
cout << "Failed testIntegrateMixedBasis with maxDiff = " << maxDiff << endl;
}
return success;
}
开发者ID:CamelliaDPG,项目名称:Camellia,代码行数:92,代码来源:LinearTermTests.cpp
示例13: main
int main(int argc, char *argv[])
{
Teuchos::GlobalMPISession mpiSession(&argc, &argv, 0);
int spaceDim = 2;
int meshWidth = 2;
bool conformingTraces = true;
int H1Order = 2, delta_k = 3;
double domainWidth = 1.0e-3;
bool diagScaling = false;
double h = domainWidth / meshWidth;
double weight = h / 4.0; // ratio of area of square with sidelength h to its perimeter
double sigma_weight = 1.0; // h / 4.0; // sigma = sigma_weight * u->grad()
Space uHatSpace = conformingTraces ? HGRAD : L2;
VarFactoryPtr vf = VarFactory::varFactory();
// fields
VarPtr u = vf->fieldVar("u");
VarPtr sigma = vf->fieldVar("sigma", VECTOR_L2);
// traces
VarPtr u_hat = vf->traceVar("u_hat", uHatSpace);
VarPtr sigma_n = vf->fluxVar("sigma_n");
// tests
VarPtr v = vf->testVar("v", HGRAD);
VarPtr tau = vf->testVar("tau", HDIV);
BFPtr bf = BF::bf(vf);
// standard BF:
// bf->addTerm(sigma, v->grad());
// bf->addTerm(sigma_n, v);
//
// bf->addTerm(sigma, tau);
// bf->addTerm(u, tau->div());
// bf->addTerm(-u_hat, tau->dot_normal());
// weighted BF:
bf->addTerm(sigma, v->grad());
bf->addTerm(weight * sigma_n, v);
bf->addTerm(sigma, tau);
bf->addTerm(sigma_weight * u, tau->div());
bf->addTerm(- sigma_weight * weight * u_hat, tau->dot_normal());
IPPtr ip = IP::ip();
// standard IP:
ip->addTerm(tau + v->grad());
ip->addTerm(tau->div());
ip->addTerm(v);
ip->addTerm(tau);
// weighted IP:
// ip->addTerm(tau + v->grad());
// ip->addTerm(sigma_weight * tau->div());
// ip->addTerm(max(sigma_weight,1e-3) * v);
// ip->addTerm(sigma_weight * weight * tau);
BCPtr bc = BC::bc();
bc->addDirichlet(u_hat, SpatialFilter::allSpace(), Function::zero());
RHSPtr rhs = RHS::rhs();
rhs->addTerm(1.0 * sigma_weight * v);
vector<double> dimensions(spaceDim,domainWidth);
vector<int> elementCounts(spaceDim,meshWidth);
MeshPtr mesh = MeshFactory::rectilinearMesh(bf, dimensions, elementCounts, H1Order, delta_k);
SolutionPtr soln = Solution::solution(mesh, bc, rhs, ip);
soln->setUseCondensedSolve(true);
soln->initializeLHSVector();
soln->initializeStiffnessAndLoad();
soln->populateStiffnessAndLoad();
Teuchos::RCP<Epetra_RowMatrix> stiffness = soln->getStiffnessMatrix();
double condNumber = conditionNumberLAPACK(*stiffness, diagScaling);
cout << "condest (1-norm): " << condNumber << endl;
return 0;
}
开发者ID:vijaysm,项目名称:Camellia,代码行数:88,代码来源:ConditionNumberMinimizationExperiment.cpp
示例14: testPoisson
bool VectorizedBasisTestSuite::testPoisson()
{
bool success = true;
//////////////////// DECLARE VARIABLES ///////////////////////
// define test variables
VarFactoryPtr varFactory = VarFactory::varFactory();
VarPtr tau = varFactory->testVar("\\tau", HDIV);
VarPtr v = varFactory->testVar("v", HGRAD);
// define trial variables
VarPtr uhat = varFactory->traceVar("\\widehat{u}");
VarPtr sigma_n = varFactory->fluxVar("\\widehat{\\sigma_{n}}");
VarPtr u = varFactory->fieldVar("u");
VarPtr sigma = varFactory->fieldVar("\\sigma", VECTOR_L2);
//////////////////// DEFINE BILINEAR FORM ///////////////////////
BFPtr bf = Teuchos::rcp( new BF(varFactory) );
// tau terms:
bf->addTerm(sigma, tau);
bf->addTerm(u, tau->div());
bf->addTerm(-uhat, tau->dot_normal());
// v terms:
bf->addTerm( sigma, v->grad() );
bf->addTerm( -sigma_n, v);
//////////////////// DEFINE INNER PRODUCT(S) ///////////////////////
IPPtr ip = bf->graphNorm();
//////////////////// SPECIFY RHS ///////////////////////
RHSPtr rhs = RHS::rhs();
FunctionPtr f = Function::constant(1.0);
rhs->addTerm( f * v );
//////////////////// CREATE BCs ///////////////////////
BCPtr bc = BC::bc();
SpatialFilterPtr boundary = SpatialFilter::allSpace();
FunctionPtr zero = Function::zero();
bc->addDirichlet(uhat, boundary, zero);
//////////////////// BUILD MESH ///////////////////////
int H1Order = 3, pToAdd = 2;
// define nodes for mesh
FieldContainer<double> meshBoundary(4,2);
meshBoundary(0,0) = 0.0; // x1
meshBoundary(0,1) = 0.0; // y1
meshBoundary(1,0) = 1.0;
meshBoundary(1,1) = 0.0;
meshBoundary(2,0) = 1.0;
meshBoundary(2,1) = 1.0;
meshBoundary(3,0) = 0.0;
meshBoundary(3,1) = 1.0;
int horizontalCells = 1, verticalCells = 1;
// create a pointer to a new mesh:
Teuchos::RCP<Mesh> mesh = MeshFactory::buildQuadMesh(meshBoundary, horizontalCells, verticalCells,
bf, H1Order, H1Order+pToAdd, false);
//////////////////// SOLVE & REFINE ///////////////////////
Teuchos::RCP<Solution> solution = Teuchos::rcp( new Solution(mesh, bc, rhs, ip) );
double energyThreshold = 0.2; // for mesh refinements
RefinementStrategy refinementStrategy( solution, energyThreshold );
#ifdef USE_VTK
VTKExporter exporter(solution, mesh, varFactory);
#endif
for (int refIndex=0; refIndex<=4; refIndex++)
{
solution->solve(false);
#ifdef USE_VTK
// output commented out because it's not properly part of the test.
// stringstream outfile;
// outfile << "test_" << refIndex;
// exporter.exportSolution(outfile.str());
#endif
if (refIndex < 4)
refinementStrategy.refine(false); // don't print to console
}
return success;
}
开发者ID:CamelliaDPG,项目名称:Camellia,代码行数:84,代码来源:VectorizedBasisTestSuite.cpp
示例15: testRieszInversion
// tests Riesz inversion by integration by parts
bool LinearTermTests::testRieszInversion()
{
bool success = true;
//////////////////// DECLARE VARIABLES ///////////////////////
// define test variables
VarFactoryPtr varFactory = VarFactory::varFactory();
VarPtr tau = varFactory->testVar("\\tau", HDIV);
VarPtr v = varFactory->testVar("v", HGRAD);
// define trial variables
VarPtr uhat = varFactory->traceVar("\\widehat{u}");
VarPtr beta_n_u_minus_sigma_n = varFactory->fluxVar("\\widehat{\\beta \\cdot n u - \\sigma_{n}}");
VarPtr u = varFactory->fieldVar("u");
VarPtr sigma1 = varFactory->fieldVar("\\sigma_1");
VarPtr sigma2 = varFactory->fieldVar("\\sigma_2");
vector<double> beta;
beta.push_back(1.0);
beta.push_back(0.0);
double eps = .01;
//////////////////// DEFINE BILINEAR FORM ///////////////////////
BFPtr confusionBF = Teuchos::rcp( new BF(varFactory) );
// tau terms:
confusionBF->addTerm(sigma1 / eps, tau->x());
confusionBF->addTerm(sigma2 / eps, tau->y());
confusionBF->addTerm(u, tau->div());
confusionBF->addTerm(uhat, -tau->dot_normal());
// v terms:
confusionBF->addTerm( sigma1, v->dx() );
confusionBF->addTerm( sigma2, v->dy() );
confusionBF->addTerm( -u, beta * v->grad() );
confusionBF->addTerm( beta_n_u_minus_sigma_n, v);
//////////////////// BUILD MESH ///////////////////////
// define nodes for mesh
int H1Order = 1;
int pToAdd = 1;
FieldContainer<double> quadPoints(4,2);
quadPoints(0,0) = 0.0; // x1
quadPoints(0,1) = 0.0; // y1
quadPoints(1,0) = 1.0;
quadPoints(1,1) = 0.0;
quadPoints(2,0) = 1.0;
quadPoints(2,1) = 1.0;
quadPoints(3,0) = 0.0;
quadPoints(3,1) = 1.0;
int nCells = 1;
int horizontalCells = nCells, verticalCells = nCells;
// create a pointer to a new mesh:
Teuchos::RCP<Mesh> myMesh = MeshFactory::buildQuadMesh(quadPoints, horizontalCells, verticalCells,
confusionBF, H1Order, H1Order+pToAdd);
ElementTypePtr elemType = myMesh->getElement(0)->elementType();
BasisCachePtr basisCache = Teuchos::rcp(new BasisCache(elemType, myMesh));
vector<GlobalIndexType> cellIDs;
vector<ElementPtr> elems = myMesh->activeElements();
vector<ElementPtr>::iterator elemIt;
for (elemIt=elems.begin(); elemIt!=elems.end(); elemIt++)
{
int cellID = (*elemIt)->cellID();
cellIDs.push_back(cellID);
}
bool createSideCacheToo = true;
basisCache->setPhysicalCellNodes(myMesh->physicalCellNodesGlobal(elemType), cellIDs, createSideCacheToo);
LinearTermPtr integrand = Teuchos::rcp(new LinearTerm);// residual
LinearTermPtr integrandIBP = Teuchos::rcp(new LinearTerm);// residual
vector<double> e1(2); // (1,0)
vector<double> e2(2); // (0,1)
e1[0] = 1;
e2[1] = 1;
FunctionPtr n = Function::normal();
FunctionPtr X = Function::xn(1);
FunctionPtr Y = Function::yn(1);
FunctionPtr testFxn1 = X;
FunctionPtr testFxn2 = Y;
FunctionPtr divTestFxn = testFxn1->dx() + testFxn2->dy();
FunctionPtr vectorTest = testFxn1*e1 + testFxn2*e2;
integrand->addTerm(divTestFxn*v);
integrandIBP->addTerm(vectorTest*n*v - vectorTest*v->grad()); // boundary term
IPPtr sobolevIP = Teuchos::rcp(new IP);
sobolevIP->addTerm(v);
sobolevIP->addTerm(tau);
Teuchos::RCP<RieszRep> riesz = Teuchos::rcp(new RieszRep(myMesh, sobolevIP, integrand));
// riesz->setPrintOption(true);
riesz->computeRieszRep();
//.........这里部分代码省略.........
开发者ID:CamelliaDPG,项目名称:Camellia,代码行数:101,代码来源:LinearTermTests.cpp
示例16: main
int main(int argc, char *argv[])
{
#ifdef HAVE_MPI
Teuchos::GlobalMPISession mpiSession(&argc, &argv,0);
Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
Epetra_SerialComm Comm;
#endif
int commRank = Teuchos::GlobalMPISession::getRank();
Comm.Barrier(); // set breakpoint here to allow debugger attachment to other MPI processes than the one you automatically attached to.
Teuchos::CommandLineProcessor cmdp(false,true); // false: don't throw exceptions; true: do return errors for unrecognized options
// problem parameters:
double mu = 0.1;
double permCoef = 1e4;
int numRefs = 0;
int k = 2, delta_k = 2;
string norm = "Graph";
cmdp.setOption("polyOrder",&k,"polynomial order for field variable u");
cmdp.setOption("delta_k", &delta_k, "test space polynomial order enrichment");
cmdp.setOption("numRefs",&numRefs,"number of refinements");
cmdp.setOption("norm", &norm, "norm");
cmdp.setOption("mu", &mu, "mu");
cmdp.setOption("permCoef", &permCoef, "Permeability coefficient");
if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL)
{
#ifdef HAVE_MPI
MPI_Finalize();
#endif
return -1;
}
FunctionPtr zero = TFunction<double>::zero();
FunctionPtr one = TFunction<double>::constant(1);
FunctionPtr sin2pix = Teuchos::rcp( new Sin_ax(2*pi) );
FunctionPtr cos2pix = Teuchos::rcp( new Cos_ax(2*pi) );
FunctionPtr sin2piy = Teuchos::rcp( new Sin_ay(2*pi) );
FunctionPtr cos2piy = Teuchos::rcp( new Cos_ay(2*pi) );
FunctionPtr u1_exact = sin2pix*cos2piy;
FunctionPtr u2_exact = -cos2pix*sin2piy;
FunctionPtr x2 = TFunction<double>::xn(2);
FunctionPtr y2 = TFunction<double>::yn(2);
FunctionPtr p_exact = x2*y2 - 1./9;
FunctionPtr permInv = permCoef*(sin2pix + 1.1);
VarFactoryPtr vf = VarFactory::varFactory();
//fields:
VarPtr sigma1 = vf->fieldVar("sigma1", VECTOR_L2);
VarPtr sigma2 = vf->fieldVar("sigma2", VECTOR_L2);
VarPtr u1 = vf->fieldVar("u1", L2);
VarPtr u2 = vf->fieldVar("u2", L2);
VarPtr p = vf->fieldVar("p", L2);
// traces:
VarPtr u1hat = vf->traceVar("u1hat");
VarPtr u2hat = vf->traceVar("u2hat");
VarPtr t1c = vf->fluxVar("t1c");
VarPtr t2c = vf->fluxVar("t2c");
// test:
VarPtr v1 = vf->testVar("v1", HGRAD);
VarPtr v2 = vf->testVar("v2", HGRAD);
VarPtr tau1 = vf->testVar("tau1", HDIV);
VarPtr tau2 = vf->testVar("tau2", HDIV);
VarPtr q = vf->testVar("q", HGRAD);
BFPtr bf = Teuchos::rcp( new BF(vf) );
bf->addTerm(1./mu*sigma1, tau1);
bf->addTerm(1./mu*sigma2, tau2);
bf->addTerm(u1, tau1->div());
bf->addTerm(u2, tau2->div());
bf->addTerm(-u1hat, tau1->dot_normal());
bf->addTerm(-u2hat, tau2->dot_normal());
bf->addTerm(sigma1, v1->grad());
bf->addTerm(sigma2, v2->grad());
bf->addTerm(-p, v1->dx());
bf->addTerm(-p, v2->dy());
bf->addTerm(t1c, v1);
bf->addTerm(t2c, v2);
bf->addTerm(mu*permInv*u1, v1);
bf->addTerm(mu*permInv*u2, v2);
bf->addTerm(-u1, q->dx());
bf->addTerm(-u2, q->dy());
bf->addTerm(u1hat, q->times_normal_x());
bf->addTerm(u2hat, q->times_normal_y());
RHSPtr rhs = RHS::rhs();
BCPtr bc = BC::bc();
SpatialFilterPtr y_equals_one = SpatialFilter::matchingY(1.0);
//.........这里部分代码省略.........
开发者ID:CamelliaDPG,项目名称:Camellia,代码行数:101,代码来源:BrinkmanDriver.cpp
示例17: phi_hat
PoissonFormulation::PoissonFormulation(int spaceDim, bool useConformingTraces, PoissonFormulationChoice formulationChoice)
{
_spaceDim = spaceDim;
if (formulationChoice == ULTRAWEAK)
{
Space tauSpace = (spaceDim > 1) ? HDIV : HGRAD;
Space phi_hat_space = useConformingTraces ? HGRAD : L2;
Space psiSpace = (spaceDim > 1) ? VECTOR_L2 : L2;
// fields
VarPtr phi;
VarPtr psi;
// traces
VarPtr phi_hat, psi_n_hat;
// tests
VarPtr q;
VarPtr tau;
VarFactoryPtr vf = VarFactory::varFactory();
phi = vf->fieldVar(S_PHI);
psi = vf->fieldVar(S_PSI, psiSpace);
TFunctionPtr<double> parity = TFunction<double>::sideParity();
if (spaceDim > 1)
phi_hat = vf->traceVar(S_PHI_HAT, phi, phi_hat_space);
else
phi_hat = vf->fluxVar(S_PHI_HAT, phi * (Function::normal_1D() * parity), phi_hat_space); // for spaceDim==1, the "normal" component is in the flux-ness of phi_hat (it's a plus or minus 1)
TFunctionPtr<double> n = TFunction<double>::normal();
if (spaceDim > 1)
psi_n_hat = vf->fluxVar(S_PSI_N_HAT, psi * (n * parity));
else
psi_n_hat = vf->fluxVar(S_PSI_N_HAT, psi * (Function::normal_1D() * parity));
q = vf->testVar(S_Q, HGRAD);
tau = vf->testVar(S_TAU, tauSpace);
_poissonBF = Teuchos::rcp( new BF(vf) );
if (spaceDim==1)
{
// for spaceDim==1, the "normal" component is in the flux-ness of phi_hat (it's a plus or minus 1)
_poissonBF->addTerm(phi, tau->dx());
_poissonBF->addTerm(psi, tau);
_poissonBF->addTerm(-phi_hat, tau);
_poissonBF->addTerm(-psi, q->dx());
_poissonBF->addTerm(psi_n_hat, q);
}
else
{
_poissonBF->addTerm(phi, tau->div());
_poissonBF->addTerm(psi, tau);
_poissonBF->addTerm(-phi_hat, tau->dot_normal());
_poissonBF->addTerm(-psi, q->grad());
_poissonBF->addTerm(psi_n_hat, q);
}
}
else if ((formulationChoice == PRIMAL) || (formulationChoice == CONTINUOUS_GALERKIN))
{
// field
VarPtr phi;
// flux
VarPtr psi_n_hat;
// tests
VarPtr q;
VarFactoryPtr vf = VarFactory::varFactory();
phi = vf->fieldVar(S_PHI, HGRAD);
TFunctionPtr<double> parity = TFunction<double>::sideParity();
TFunctionPtr<double> n = TFunction<double>::normal();
if (formulationChoice == PRIMAL)
{
if (spaceDim > 1)
psi_n_hat = vf->fluxVar(S_PSI_N_HAT, phi->grad() * (n * parity));
else
psi_n_hat = vf->fluxVar(S_PSI_N_HAT, phi->dx() * (Function::normal_1D() * parity));
}
q = vf->testVar(S_Q, HGRAD);
_poissonBF = BF::bf(vf);
_poissonBF->addTerm(-phi->grad(), q->grad());
if (formulationChoice == CONTINUOUS_GALERKIN)
{
FunctionPtr boundaryIndicator = Function::meshBoundaryCharacteristic();
_poissonBF->addTerm(phi->grad() * n, boundaryIndicator * q);
}
else // primal
{
//.........这里部分代码省略.........
开发者ID:CamelliaDPG,项目名称:Camellia,代码行数:101,代码来源:PoissonFormulation.cpp
示例18: main
//.........这里部分代码省略.........
// Domain from Evans Hughes for Navier-Stokes:
// quadPointsKovasznay(0,0) = 0.0; // x1
// quadPointsKovasznay(0,1) = -0.5; // y1
// quadPointsKovasznay(1,0) = 1.0;
// quadPointsKovasznay(1,1) = -0.5;
// quadPointsKovasznay(2,0) = 1.0;
// quadPointsKovasznay(2,1) = 0.5;
// quadPointsKovasznay(3,0) = 0.0;
// quadPointsKovasznay(3,1) = 0.5;
// double Re = 10.0; // Cockburn Kanschat Stokes
// double Re = 40.0; // Evans Hughes Navier-Stokes
// double Re = 1000.0;
string formulationTypeStr = "vgp";
FunctionPtr u1_exact, u2_exact, p_exact;
int numCellsFineMesh = 20; // for computing a zero-mean pressure
int H1OrderFineMesh = 5;
// VGPNavierStokesProblem(double Re, Intrepid::FieldContainer<double> &quadPoints, int horizontalCells,
// int verticalCells, int H1Order, int pToAdd,
// TFunctionPtr<double> u1_0, TFunctionPtr<double> u2_0, TFunctionPtr<double> f1, TFunctionPtr<double> f2,
// bool enrichVelocity = false, bool enhanceFluxes = false)
FunctionPtr zero = Function::zero();
VGPNavierStokesProblem zeroProblem = VGPNavierStokesProblem(Re, quadPointsKovasznay,
numCellsFineMesh, numCellsFineMesh,
H1OrderFineMesh, pToAdd,
zero, zero, zero, useCompliantNorm || useStokesCompliantNorm);
VarFactoryPtr varFactory = VGPStokesFormulation::vgpVarFactory();
VarPtr u1_vgp = varFactory->fieldVar(VGP_U1_S);
VarPtr u2_vgp = varFactory->fieldVar(VGP_U2_S);
VarPtr sigma11_vgp = varFactory->fieldVar(VGP_SIGMA11_S);
VarPtr sigma12_vgp = varFactory->fieldVar(VGP_SIGMA12_S);
VarPtr sigma21_vgp = varFactory->fieldVar(VGP_SIGMA21_S);
VarPtr sigma22_vgp = varFactory->fieldVar(VGP_SIGMA22_S);
VarPtr p_vgp = varFactory->fieldVar(VGP_P_S);
VarPtr v1_vgp = varFactory->testVar(VGP_V1_S, HGRAD);
VarPtr v2_vgp = varFactory->testVar(VGP_V2_S, HGRAD);
// if (rank==0) {
// cout << "bilinear form with zero background flow:\n";
// zeroProblem.bf()->printTrialTestInteractions();
// }
VGPStokesFormulation stokesForm(1/Re);
NavierStokesFormulation::setKovasznay(Re, zeroProblem.mesh(), u1_exact, u2_exact, p_exact);
// if (rank==0) cout << "Stokes bilinearForm: " << stokesForm.bf()->displayString() << endl;
map< string, string > convergenceDataForMATLAB; // key: field file name
for (int polyOrder = minPolyOrder; polyOrder <= maxPolyOrder; polyOrder++)
{
int H1Order = polyOrder + 1;
int numCells1D = pow(2.0,minLogElements);
if (rank==0)
开发者ID:CamelliaDPG,项目名称:Camellia,代码行数:67,代码来源:NavierStokesStudy.cpp
示例19: setup
void FunctionTests::setup()
{
//////////////////// DECLARE VARIABLES ///////////////////////
// define test variables
VarFactoryPtr varFactory = VarFactory::varFactory();
VarPtr tau = varFactory->testVar("\\tau", HDIV);
VarPtr v = varFactory->testVar("v", HGRAD);
// define trial variables
VarPtr uhat = varFactory->traceVar("\\widehat{u}");
VarPtr beta_n_u_minus_sigma_n = varFactory->fluxVar("\\widehat{\\beta \\cdot n u - \\sigma_{n}}");
VarPtr u = varFactory->fieldVar("u");
VarPtr sigma1 = varFactory->fieldVar("\\sigma_1");
VarPtr sigma2 = varFactory->fieldVar("\\sigma_2");
vector<double> beta_const;
beta_const.push_back(2.0);
beta_const.pu
|
请发表评论