Skip to content
Snippets Groups Projects
Commit e11f58ce authored by Müller, Felix's avatar Müller, Felix
Browse files

Merge branch 'feature/lagrangebasis_with_single_node' into 'master'

Add overload for LagrangeBasis with single node

See merge request !152
parents eda87e6c d9715c4a
No related branches found
No related tags found
1 merge request!152Add overload for LagrangeBasis with single node
Showing with 73 additions and 71 deletions
......@@ -35,4 +35,4 @@ endif ()
if (dune-foamgrid_FOUND)
add_amdis_executable(NAME surface.2d SOURCES surface.cc DIM 2 DOW 3 ALBERTA_GRID)
add_dependencies(examples surface.2d)
endif ()
\ No newline at end of file
endif ()
......@@ -11,10 +11,8 @@
#include <amdis/Integrate.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/LocalOperators.hpp>
#include <amdis/common/Literals.hpp>
using namespace AMDiS;
using namespace Dune::Indices;
// 1 component with polynomial degree 2
using Param = YaspGridBasis<GRIDDIM, 2>;
......@@ -28,7 +26,7 @@ void run(SetBoundary setBoundary)
setBoundary(prob.boundaryManager());
auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
prob.addMatrixOperator(opL, 0, 0);
prob.addMatrixOperator(opL);
auto f = [](auto const& x) {
double r2 = dot(x,x);
......@@ -36,18 +34,18 @@ void run(SetBoundary setBoundary)
return -(400.0 * r2 - 20.0 * x.size()) * ux;
};
auto opForce = makeOperator(tag::test{}, f, 6);
prob.addVectorOperator(opForce, 0);
prob.addVectorOperator(opForce);
// set boundary condition
auto g = [](auto const& x){ return std::exp(-10.0 * dot(x,x)); };
prob.addDirichletBC(BoundaryType{1}, 0, 0, g);
prob.addDirichletBC(BoundaryType{1}, g);
AdaptInfo adaptInfo("adapt");
prob.assemble(adaptInfo);
prob.solve(adaptInfo);
double errorL2 = integrate(sqr(g - prob.solution(0)), prob.gridView(), 6);
double errorL2 = integrate(sqr(g - prob.solution()), prob.gridView(), 6);
msg("error_L2 = {}", errorL2);
}
......@@ -69,17 +67,17 @@ void run_periodic()
prob.boundaryManager().setBoxBoundary({-1,-1,1,1});
auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
prob.addMatrixOperator(opL, 0, 0);
prob.addMatrixOperator(opL);
auto opForce = makeOperator(tag::test{}, 1.0, 6);
prob.addVectorOperator(opForce, 0);
prob.addVectorOperator(opForce);
// set boundary condition
auto g = [](auto const& x) {
return std::sin(0.5*M_PI*x[1])*std::sin(2*M_PI * (0.25 + x[0]))
+ std::cos(0.5*M_PI*x[1])*std::sin(2*M_PI * (-0.25 + x[0]));
};
prob.addDirichletBC(BoundaryType{1}, 0, 0, g);
prob.addDirichletBC(BoundaryType{1}, g);
typename ElliptProblem::WorldMatrix A{{1.0,0.0}, {0.0,1.0}};
typename ElliptProblem::WorldVector b{1.0, 0.0};
......
......@@ -6,7 +6,6 @@
#include <amdis/AMDiS.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/localoperators/ConvectionDiffusionOperator.hpp>
#include <amdis/common/Literals.hpp>
using namespace AMDiS;
......@@ -19,20 +18,18 @@ int main(int argc, char** argv)
{
Environment env(argc, argv);
using namespace Dune::Indices;
ElliptProblem prob("ellipt");
prob.initialize(INIT_ALL);
// -div(A*grad(u)) + div(b*u) + c*u = f
auto opCD = convectionDiffusion(/*A=*/1.0, /*b=*/0.0, /*c=*/1.0, /*f=*/1.0);
prob.addMatrixOperator(opCD, 0, 0);
prob.addVectorOperator(opCD, 0);
prob.addMatrixOperator(opCD);
prob.addVectorOperator(opCD);
// set boundary condition
auto predicate = [](auto const& x){ return x[0] < 1.e-8 || x[1] < 1.e-8; }; // define boundary
auto dbcValues = [](auto const& x){ return 0.0; }; // set value
prob.addDirichletBC(predicate, 0, 0, dbcValues);
prob.addDirichletBC(predicate, dbcValues);
AdaptInfo adaptInfo("adapt");
prob.assemble(adaptInfo);
......
......@@ -8,7 +8,6 @@
using Grid = Dune::YaspGrid<GRIDDIM>;
using namespace AMDiS;
using namespace Dune::Indices;
int main(int argc, char** argv)
{
......@@ -34,13 +33,13 @@ int main(int argc, char** argv)
prob.initialize(INIT_ALL);
auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
prob.addMatrixOperator(opL, 0, 0);
prob.addMatrixOperator(opL);
auto opForce = makeOperator(tag::test{}, f, 6);
prob.addVectorOperator(opForce, 0);
prob.addVectorOperator(opForce);
// set boundary condition
prob.addDirichletBC(1, 0, 0, g);
prob.addDirichletBC(1, g);
AdaptInfo adaptInfo("adapt");
......@@ -65,9 +64,9 @@ int main(int argc, char** argv)
prob.assemble(adaptInfo);
prob.solve(adaptInfo);
double errorL2 = integrate(sqr(g - prob.solution(0)), prob.gridView(), 6);
double errorL2 = integrate(sqr(g - prob.solution()), prob.gridView(), 6);
errL2.push_back(std::sqrt(errorL2));
double errorH1 = errorL2 + integrate(unary_dot(grad_g - gradientAtQP(prob.solution(0))), prob.gridView(), 6);
double errorH1 = errorL2 + integrate(unary_dot(grad_g - gradientAtQP(prob.solution())), prob.gridView(), 6);
errH1.push_back(std::sqrt(errorH1));
}
......
......@@ -9,7 +9,6 @@
#include <amdis/ProblemInstat.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/GridFunctions.hpp>
#include <amdis/common/Literals.hpp>
using namespace AMDiS;
......@@ -34,23 +33,23 @@ int main(int argc, char** argv)
auto invTau = std::ref(probInstat.invTau());
auto opTimeLhs = makeOperator(tag::test_trial{}, invTau);
prob.addMatrixOperator(opTimeLhs, 0, 0);
prob.addMatrixOperator(opTimeLhs);
auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
prob.addMatrixOperator(opL, 0, 0);
prob.addMatrixOperator(opL);
auto opTimeRhs = makeOperator(tag::test{},
invokeAtQP([invTau](double u) { return u * invTau.get(); }, prob.solution(0)), 2);
prob.addVectorOperator(opTimeRhs, 0);
invokeAtQP([invTau](double u) { return u * invTau.get(); }, prob.solution()), 2);
prob.addVectorOperator(opTimeRhs);
auto opForce = makeOperator(tag::test{}, [](auto const& x) { return -1.0; }, 0);
prob.addVectorOperator(opForce, 0);
prob.addVectorOperator(opForce);
// set boundary condition
auto predicate = [](auto const& p){ return p[0] < 1.e-8 || p[1] < 1.e-8; };
auto dbcValues = [](auto const& p){ return 0.0; };
prob.addDirichletBC(predicate, 0, 0, dbcValues);
prob.addDirichletBC(predicate, dbcValues);
AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo);
adapt.adapt();
......
......@@ -8,7 +8,7 @@ ellipt->solver->relative tolerance: 1.e-10
ellipt->solver->info: 1
ellipt->solver->precon: ilu
ellipt->output[0]->format: vtk
ellipt->output[0]->filename: boundary.2d
ellipt->output[0]->name: u
ellipt->output[0]->output directory: ./output
ellipt->output->format: vtk
ellipt->output->filename: boundary.2d
ellipt->output->name: u
ellipt->output->output directory: ./output
......@@ -16,7 +16,7 @@ ch->solver: direct
ch->solver->max iteration: 1000
ch->solver->relative tolerance: 1e-6
ch->solver->break if tolerance not reached: 1
ch->solver->info: 2
ch->solver->info: 0
ch->output[0]->format: vtk
ch->output[0]->filename: phi.2d
......@@ -34,4 +34,4 @@ ch->output[1]->animation: 1
adapt->timestep: 0.001
adapt->start time: 0.0
adapt->end time: 1.0
adapt->end time: 0.5
......@@ -13,7 +13,7 @@ ellipt->solver->max iteration: 10000
ellipt->solver->relative tolerance: 1.e-10
ellipt->solver->precon: ilu
ellipt->output[0]->format: vtk backup
ellipt->output[0]->filename: ellipt.2d
ellipt->output[0]->name: u
ellipt->output[0]->output directory: ./output
ellipt->output->format: vtk backup
ellipt->output->filename: ellipt.2d
ellipt->output->name: u
ellipt->output->output directory: ./output
......@@ -13,7 +13,7 @@ ellipt->solver->max iteration: 10000
ellipt->solver->relative tolerance: 1.e-10
ellipt->solver->precon: ilu
ellipt->output[0]->format: vtk
ellipt->output[0]->filename: ellipt.3d
ellipt->output[0]->name: u
ellipt->output[0]->output directory: ./output
ellipt->output->format: vtk
ellipt->output->filename: ellipt.3d
ellipt->output->name: u
ellipt->output->output directory: ./output
......@@ -12,12 +12,12 @@ heat->solver->absolute tolerance: 1e-6
heat->solver->info: 1
heat->solver->precon: ilu
heat->output[0]->format: vtk backup
heat->output[0]->name: u
heat->output[0]->filename: heat.2d
heat->output[0]->output directory: output
heat->output[0]->mode: 1
heat->output[0]->animation: 1
heat->output->format: vtk backup
heat->output->name: u
heat->output->filename: heat.2d
heat->output->output directory: output
heat->output->mode: 1
heat->output->animation: 1
adapt->timestep: 0.1
adapt->start time: 0.0
......
......@@ -9,7 +9,7 @@ surface->solver->relative tolerance: 1.e-6
surface->solver->info: 1
surface->solver->precon: ilu
surface->output[0]->format: vtk
surface->output[0]->filename: surface.2d
surface->output[0]->name: u
surface->output[0]->output directory: ./output
surface->output->format: vtk
surface->output->filename: surface.2d
surface->output->name: u
surface->output->output directory: ./output
......@@ -15,10 +15,8 @@
#include <amdis/Integrate.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/LocalOperators.hpp>
#include <amdis/common/Literals.hpp>
using namespace AMDiS;
using namespace Dune::Indices;
template <class Grid>
void run(Grid& grid)
......@@ -28,18 +26,18 @@ void run(Grid& grid)
using Traits = LagrangeBasis<Grid, 2>;
ProblemStat<Traits> prob("ellipt", grid);
prob.initialize(INIT_ALL);
prob.boundaryManager().setBoxBoundary({1,2,2,1});
prob.boundaryManager()->setBoxBoundary({1,2,2,1});
auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
prob.addMatrixOperator(opL, 0, 0);
prob.addMatrixOperator(opL);
// set dirichlet boundary condition
auto g = [](auto const& x) { return 0.0; };
prob.addDirichletBC(BoundaryType{1}, 0, 0, g);
prob.addDirichletBC(BoundaryType{1}, g);
// set neumann boundary condition
auto opNeumann = makeOperator(tag::test{}, 1.0);
prob.addVectorOperator(BoundaryType{2}, opNeumann, 0);
prob.addVectorOperator(BoundaryType{2}, opNeumann);
AdaptInfo adaptInfo("adapt");
......
......@@ -11,7 +11,6 @@
#include <amdis/Integrate.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/LocalOperators.hpp>
#include <amdis/common/Literals.hpp>
using namespace AMDiS;
......
......@@ -3,7 +3,6 @@
#include <amdis/Integrate.hpp>
#include <amdis/LocalOperators.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/common/Literals.hpp>
#include <dune/grid/albertagrid/albertareader.hh>
#include <dune/foamgrid/foamgrid.hh>
......@@ -12,7 +11,6 @@
#include "SphereMapping.hpp"
using namespace AMDiS;
using namespace Dune::Indices;
using Grid = Dune::FoamGrid<2,3>;
using Basis = LagrangeBasis<Grid, 1>;
......@@ -44,14 +42,14 @@ int main(int argc, char** argv)
prob.initialize(INIT_ALL);
auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
prob.addMatrixOperator(opL, 0, 0);
prob.addMatrixOperator(opL);
double kappa = Parameters::get<double>("surface->kappa").value_or(1.0);
auto opM = makeOperator(tag::test_trial{}, -kappa*kappa);
prob.addMatrixOperator(opM, 0, 0);
prob.addMatrixOperator(opM);
auto opForce = makeOperator(tag::test{}, X(0) + X(1) + X(2));
prob.addVectorOperator(opForce, 0);
prob.addVectorOperator(opForce);
AdaptInfo adaptInfo("adapt");
prob.assemble(adaptInfo);
......
......@@ -266,6 +266,12 @@ namespace AMDiS
RowTreePath row, ColTreePath col,
Values const& values);
template <class Identifier, class Values>
void addDirichletBC(Identifier&& id, Values&& values)
{
addDirichletBC(FWD(id), RootTreePath{}, RootTreePath{}, FWD(values));
}
/// Add a periodic boundary conditions to the system, by specifying a face transformation
/// y = A*x + b of coordinates. We assume, that A is orthonormal.
void addPeriodicBC(BoundaryType id, WorldMatrix const& A, WorldVector const& b);
......
......@@ -18,6 +18,17 @@ namespace AMDiS
template <bool same, int... degs>
struct LagrangePreBasisCreatorImpl;
// specialization for single node basis
template <int deg>
struct LagrangePreBasisCreatorImpl<true, deg>
{
static auto create()
{
using namespace Dune::Functions::BasisFactory;
return lagrange<deg>();
}
};
// specialization for composite basis
template <int... degs>
struct LagrangePreBasisCreatorImpl<false, degs...>
......
......@@ -7,7 +7,6 @@
#include <amdis/Integrate.hpp>
#include <amdis/LocalOperators.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/common/Literals.hpp>
#include <amdis/common/FieldMatVec.hpp>
......@@ -20,12 +19,10 @@ int main(int argc, char** argv)
{
Environment env(argc, argv);
using namespace Dune::Indices;
ElliptProblem prob("ellipt");
prob.initialize(INIT_ALL);
auto u = prob.solution(_0);
auto u = prob.solution();
// eval a functor at global coordinates (at quadrature points)
auto expr1 = evalAtQP([](Dune::FieldVector<double, 2> const& x) { return x[0] + x[1]; });
......@@ -42,7 +39,7 @@ int main(int argc, char** argv)
auto expr8 = X(0);
// Evaluation of the DOFVector (component) at quadrature points
auto expr9 = prob.solution(_0);
auto expr9 = prob.solution();
// ---------------------------------
......
......@@ -23,7 +23,7 @@ int main(int argc, char** argv)
ElliptProblem prob("ellipt");
prob.initialize(INIT_ALL);
auto u = prob.solution(0);
auto u = prob.solution();
auto gv = u.basis().gridView();
u << 1.0;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment