Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifdef HAVE_CONFIG_H
# include "config.h"
#endif
#if ! HAVE_DUNE_FUNCTIONS
#error "Require dune-functions!"
#endif
#if ! HAVE_DUNE_ALUGRID
#error "Require dune-alugrid!"
#endif
#if ! HAVE_ALBERTA
#error "Require Alberta!"
#endif
#include <iostream>
#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
#include <dune/common/exceptions.hh> // We use exceptions
#include <dune/alugrid/grid.hh>
#include <dune/multimesh/mmhierarchiciterator.hh>
#include <dune/multimesh/mmgridfactory.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/grid/albertagrid/albertareader.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh>
#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/mtl/interface/umfpack_solve.hpp>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
Praetorius, Simon
committed
#include "phasefield4.hh"
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
using namespace Dune;
using namespace Dune::Indices;
struct MTLVector
{
using value_type = double;
MTLVector(mtl::dense_vector<double>& vec)
: vec_(vec) {}
void resize(std::size_t s) { vec_.change_dim(s); }
std::size_t size() const { return mtl::size(vec_); }
decltype(auto) operator[](std::size_t i) { return vec_[i]; }
decltype(auto) operator[](std::size_t i) const { return vec_[i]; }
decltype(auto) operator[](ReservedVector<long unsigned int, 1> i) { return vec_[i[0]]; }
decltype(auto) operator[](ReservedVector<long unsigned int, 1> i) const { return vec_[i[0]]; }
mtl::dense_vector<double>& vec_;
};
template <class VectorType, class Basis>
void writeFile(VectorType& x, Basis const& basis, std::string filename)
{
auto gv = basis.gridView();
auto xWrapper = MTLVector(x);
auto xFct = Functions::makeDiscreteGlobalBasisFunction<double>(basis, TypeTree::hybridTreePath(), xWrapper);
VTKWriter<decltype(gv)> vtkWriter(gv);
vtkWriter.addVertexData(xFct, VTK::FieldInfo("u", VTK::FieldInfo::Type::scalar, 1));
vtkWriter.write(filename);
}
int main(int argc, char** argv)
{
MPIHelper::instance(argc, argv);
assert(argc > 1);
std::string filename = argv[1];
using Grid = Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>;
AlbertaReader<Grid> albertaReader;
GridFactory<Grid> factory;
albertaReader.readGrid(filename, factory);
std::unique_ptr<Grid> gridPtr(factory.createGrid());
auto& grid = *gridPtr;
double eps = 0.05;
auto signedDistFct = [](auto const& x) { return x.two_norm() - 1.0; };
grid.globalRefine(6);
// refine grid[1] along phase-field interface
for (int i = 0; i < 8; ++i) {
std::cout << "refine " << i << "...\n";
for (auto const& elem : elements(grid.leafGridView())) {
auto geo = elem.geometry();
bool refine = false;
for (int j = 0; j < geo.corners(); ++j)
refine = refine || std::abs(signedDistFct(geo.corner(j))) < std::max(1,6-i)*eps;
if (refine)
grid.mark(1, elem);
}
grid.preAdapt();
grid.adapt();
grid.postAdapt();
}
using namespace Functions::BasisBuilder;
auto fine_basis = makeBasis(grid.leafGridView(), lagrange<1>());
using VectorType = mtl::dense_vector<double>;
using MatrixType = mtl::compressed2D<double>;
// interpolate phase-field to fine grid
VectorType phase(fine_basis.dimension());
auto phaseWrapper = MTLVector(phase);
interpolate(fine_basis, phaseWrapper, [eps,d=signedDistFct](auto const& x)
{
return 0.5*(1.0 - std::tanh(3.0*d(x)/eps));
});
writeFile(phase, fine_basis, "phase");
// assemble a sequence of systems for finer and finer grids
for (int l = 2; l < 7; ++l) {
auto coarse_basis = makeBasis(grid.levelGridView(l), lagrange<1>());
VectorType rhs(coarse_basis.dimension());
MatrixType matrix(coarse_basis.dimension(), coarse_basis.dimension());
assembleSystem(coarse_basis, fine_basis, phase, matrix, rhs, eps);
VectorType u(coarse_basis.dimension());
u = 0.0;
umfpack_solve(matrix, u, rhs);
writeFile(u, coarse_basis, "u_" + std::to_string(l));
}
}