Commit 390aede1 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

alugrid example added

parent 120d8512
......@@ -22,3 +22,4 @@ error.txt
*.out.*
*.tmp
*.*-swp
*.log
CMAKE_FLAGS=" \
-DCMAKE_CXX_COMPILER=clang++-3.8 -DCMAKE_C_COMPILER=clang-3.8 \
-DCMAKE_INSTALL_PREFIX:PATH=/opt/software/dune/dune-dec \
-DCMAKE_BUILD_TYPE=Debug \
-DCMAKE_BUILD_TYPE=RelWithDebInfo \
-DEIGEN3_INCLUDE_DIR:PATH=/opt/software/eigen/include/eigen3"
......@@ -6,4 +6,4 @@ Module: dune-dec
Version: 0.1
Maintainer: simon.praetorius@tu-dresden.de
Depends: dune-common dune-geometry dune-grid dune-localfunctions dune-functions
Suggests: dune-foamgrid dune-uggrid
\ No newline at end of file
Suggests: dune-foamgrid dune-uggrid dune-alugrid
\ No newline at end of file
......@@ -5,6 +5,7 @@
#include <dune/dec/Dec.hpp>
#include <dune/dec/GridConcepts.hpp>
#include <dune/dec/common/ConceptsBase.hpp>
#include <dune/dec/common/ScalarTypes.hpp>
#include <dune/dec/common/Output.hpp>
#include <dune/dec/linear_algebra/DOFVectorBase.hpp>
......@@ -149,6 +150,14 @@ namespace Dec
return *this;
}
template <class S,
REQUIRES(concepts::Arithmetic<S>)>
DOFVector& operator=(S const& scalar)
{
Super::setConstant(scalar);
return *this;
}
// ---------------------------------------------------------------------------------------------
/// Copy/Move assignment operator
......
......@@ -11,6 +11,9 @@ target_compile_definitions(domain_decomposition PRIVATE DEC_HAS_MPI=1)
add_dec_executable(entity_owner entity_owner.cpp)
target_compile_definitions(entity_owner PRIVATE DEC_HAS_MPI=1)
add_dec_executable(alugrid ALBERTA alugrid.cpp)
target_compile_definitions(alugrid PRIVATE DEC_HAS_MPI=1)
add_dec_executable(geometry ALBERTA geometry.cpp)
add_dec_executable(laplace ALBERTA laplace.cpp)
add_dec_executable(laplace_operator ALBERTA laplace_operator.cpp)
......@@ -21,7 +24,7 @@ add_dec_executable(transfer ALBERTA transfer.cpp)
add_dec_executable(vecellipt ALBERTA vecellipt.cpp)
# add_dec_executable(weighted_triangulation ALBERTA weighted_triangulation.cpp)
add_dec_executable(weighted_triangulation2 ALBERTA weighted_triangulation2.cpp)
add_dependencies(examples geometry laplace laplace_operator simple_grid heat weighted_triangulation2 orientation transfer domain_decomposition entity_owner)
add_dependencies(examples geometry laplace laplace_operator simple_grid heat weighted_triangulation2 orientation transfer domain_decomposition entity_owner alugrid)
add_dec_executable(helmholtz ALBERTA helmholtz.cpp)
add_dec_executable(consistency ALBERTA consistency.cpp)
......
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <list>
#include <dune/dec/DecGrid.hpp>
// #include <dune/dec/PartitionTypeMapper.hpp>
#include <dune/dec/common/Output.hpp>
#include <dune/dec/operators/Laplace.hpp>
#include <dune/dec/parallel/DistributedVector.hpp>
#include <dune/dec/parallel/DistributedMatrix.hpp>
#include <dune/dec/linear_algebra/krylov/bcgs.hpp>
#include <dune/dec/linear_algebra/iteration.hpp>
#include <dune/dec/utility/Timer.hpp>
#include <dune/grid/utility/structuredgridfactory.hh>
#include <dune/grid/utility/parmetisgridpartitioner.hh>
#include <dune/alugrid/grid.hh>
#include <dune/grid/albertagrid/albertareader.hh>
#include <dune/grid/io/file/vtk.hh>
using namespace Dec;
template <int cd, class GridView, class Predicate>
std::list<std::size_t> get_boundary(GridView const& gv, Predicate pred)
{
std::list<std::size_t> indices;
for (auto const& e : entities(gv, Codim_<cd>, Dune::Partitions::InteriorBorder{}))
{
if (pred(gv.center(e)))
indices.push_back(e.index());
}
return indices;
}
template <class GridBase>
class LoadBalance
{
public:
LoadBalance(GridBase const& grid, Dune::MPIHelper& mpihelper)
: grid_(grid)
, part_(Dune::ParMetisGridPartitioner<typename GridBase::LeafGridView>::partition(grid.leafGridView(), mpihelper))
{}
template <class Entity>
int operator()(Entity const& e) const
{
auto gv = grid_.leafGridView();
return part_[gv.indexSet().index(e)];
}
bool importRanks(std::set<int>& ranks) const { return false; }
private:
GridBase const& grid_;
std::vector<unsigned> part_;
};
int main(int argc, char** argv)
{
auto& mpihelper = Dune::MPIHelper::instance(argc, argv);
assert_msg( argc > 1, "usage: ./alugrid grid-filename [nSteps]");
Timer t;
using GridBase = Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>;
#if 1
// Create ug grid from structured grid
std::array<unsigned int, 2> n = {{10, 10}};
Dune::FieldVector<double, 2> lower = {{0, 0}};
Dune::FieldVector<double, 2> upper = {{1, 1}};
std::shared_ptr<GridBase> gridBase = Dune::StructuredGridFactory<GridBase>::createSimplexGrid(lower, upper, n);
#endif
#if 0
Dune::GridFactory<GridBase> gridFactory;
if (mpihelper.rank() == 0)
Dune::AlbertaReader<GridBase>().readGrid(argv[1], gridFactory);
std::unique_ptr<GridBase> gridBase( gridFactory.createGrid() );
#endif
msg("time(create gridBase) = ",t.elapsed(),"sec");
t.reset();
LoadBalance<GridBase> lb(*gridBase, mpihelper);
gridBase->repartition(lb);
msg("#elements = ",gridBase->leafGridView().size(0));
msg("time(loadBalance) = ",t.elapsed(),"sec");
t.reset();
using Grid = DecGrid<GridBase>;
Grid grid(*gridBase);
msg("time(create DECGrid) = ",t.elapsed(),"sec");
if (argc > 2) {
t.reset();
grid.globalRefine(std::atoi(argv[2]));
msg("#elements after global refinement = ",gridBase->leafGridView().size(0));
msg("time(global refine) = ",t.elapsed(),"sec");
}
using GridView = Grid::LeafGridView;
const GridView gv = grid.leafGridView();
auto const& indexSet = gv.indexSet();
// -----------------------------------------------------------------------------------------------
#if 0
for (auto const& v : vertices(gv, Dune::Partitions::All{}) ) {
msg("v",(indexSet.owns(v)?"*":""),"(",v.index(),") ",v.partitionType(),": ",gv.center(v));
}
for (auto const& e : elements(gv, Dune::Partitions::All{}) ) {
msg_("elem",(indexSet.owns(e)?"*":""),"(",e.index(),") ",e.partitionType(),": {");
for (auto const& v : vertices(e))
msg_(v.index(),' ');
msg("}");
}
#endif
// -----------------------------------------------------------------------------------------------
t.reset();
LaplaceBeltrami<GridView> laplace(gv);
DistributedMatrix<double, GridView, 2, 2> A(gv);
laplace.build(A);
DistributedVector<double, GridView, 2> x(gv), b(gv);
b.vector() = 1.0;
x.vector() = 0.0;
if (mpihelper.rank() == 0) {
typename Grid::GlobalCoordinate x0;
for (auto const& v : vertices(gv)) {
if (v.partitionType() == Dune::InteriorEntity && indexSet.owns(v)) {
x0 = gv.center(v);
break;
}
}
auto boundary = get_boundary<2>(gv, [&x0](auto const& x) {
return two_norm(x - x0) < 1.e-10;
});
for (auto i : boundary) {
A.matrix().clear_dirichlet_row(i);
b.vector()[i] = 0.0;
}
}
msg("time(assemble) = ",t.elapsed(),"sec");
t.reset();
ResidualIteration<double> iter(5000, 1.e-7);
solver::bcgs(A,x,b, iter);
msg("time(solve) = ",t.elapsed(),"sec");
msg("|A*x-b| = ",residuum(A,x,b));
msg("|x| = ",two_norm(x));
t.reset();
Dune::VTKWriter<GridView> vtkwriter(gv);
vtkwriter.addVertexData(x.vector(), "x");
vtkwriter.write("alugrid");
msg("time(write file) = ",t.elapsed(),"sec");
msg("grid.memory = ",grid.memory()/1024,"kB");
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment