Skip to content

Commit

Permalink
Obstacle problem working. Added some results.
Browse files Browse the repository at this point in the history
  • Loading branch information
@datafl4sh
datafl4sh committed Feb 9, 2018
1 parent 271ac6d commit 7a071b4
Show file tree
Hide file tree
Showing 6 changed files with 127 additions and 85 deletions.
1 change: 0 additions & 1 deletion CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,6 @@ set (CMAKE_MODULE_PATH
include(safeguards)

set(CMAKE_CXX_FLAGS_DEBUG "-std=c++14 -g -fpermissive -fsanitize=address")
set(CMAKE_CXX_FLAGS_DEBUG_ASAN "-std=c++14 -g -fpermissive -fsanitize=address")
set(CMAKE_CXX_FLAGS_RELEASE "-std=c++14 -O3 -mavx -g -fpermissive")
set(CMAKE_CXX_FLAGS_RELEASEASSERT "-std=c++14 -O3 -mavx -g -fpermissive")

Expand Down
71 changes: 47 additions & 24 deletions apps/obstacle/obstacle.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,13 @@ using namespace Eigen;

#include "methods/hho"


template<typename T>
struct obstacle_hho_config
{
size_t degree;
size_t max_iter;
T threshold;
};

template<typename Mesh>
void run_hho_obstacle(const Mesh& msh, size_t degree)
Expand Down Expand Up @@ -68,10 +74,10 @@ void run_hho_obstacle(const Mesh& msh, size_t degree)
};

auto sol_fun = [=](const typename Mesh::point_type& pt) -> T {
auto r = sqrt( pt.x()*pt.x() + pt.y()*pt.y() );
auto r = sqrt(pt.x()*pt.x() + pt.y()*pt.y());
auto s = r*r - r0*r0;

return std::max(s*s, 0.0);
auto t = std::max(s, 0.0);
return t*t;
};

auto bcs_fun = [&](const typename Mesh::point_type& pt) -> T {
Expand All @@ -82,15 +88,20 @@ void run_hho_obstacle(const Mesh& msh, size_t degree)
return 0.0;
};

auto num_cells = msh.cells.size();
auto num_faces = msh.faces.size();
auto fbs = face_basis<Mesh,T>::size(degree);
auto num_alpha_dofs = num_cells + fbs*num_faces;

Matrix<T, Dynamic, 1> alpha = Matrix<T, Dynamic, 1>::Zero( msh.cells.size() );
Matrix<T, Dynamic, 1> beta = Matrix<T, Dynamic, 1>::Ones( msh.cells.size() );
Matrix<T, Dynamic, 1> gamma = Matrix<T, Dynamic, 1>::Zero( msh.cells.size() );
Matrix<T, Dynamic, 1> alpha = Matrix<T, Dynamic, 1>::Zero( num_alpha_dofs );
Matrix<T, Dynamic, 1> beta = Matrix<T, Dynamic, 1>::Ones( num_cells );
Matrix<T, Dynamic, 1> gamma = Matrix<T, Dynamic, 1>::Zero( num_cells );
T c = 1.0;

size_t quadrature_degree_increase = 1;

std::vector<T> expected_solution;
expected_solution.reserve( msh.cells.size() );
expected_solution.reserve( num_cells );

size_t i = 0;
for (auto& cl : msh.cells)
Expand All @@ -105,6 +116,8 @@ void run_hho_obstacle(const Mesh& msh, size_t degree)
bool has_to_iterate = true;
while ( has_to_iterate && iter < 50 )
{

/* Prepare SILO database */
std::cout << bold << "Iteration " << iter << reset << std::endl;
std::stringstream ss;
ss << "obstacle_cycle_" << iter << ".silo";
Expand All @@ -113,7 +126,9 @@ void run_hho_obstacle(const Mesh& msh, size_t degree)
silo.create( ss.str() );
silo.add_mesh(msh, "mesh");

Matrix<T, Dynamic, 1> diff = beta + c * ( alpha - gamma );
/* Compute the beta quantity (see "Bubbles enriched quadratic finite element method for
* the 3D-elliptic obstacle problem - S. Gaddam, T. Gudi", eqn. 5.1 onwards) */
Matrix<T, Dynamic, 1> diff = beta + c * ( alpha.head(num_cells) - gamma );
std::vector<bool> in_A;
in_A.resize(diff.size());
Matrix<T, Dynamic, 1> active = Matrix<T, Dynamic, 1>::Zero(diff.size());
Expand All @@ -124,21 +139,19 @@ void run_hho_obstacle(const Mesh& msh, size_t degree)
active(i) = (diff(i) < 0);
}

/* Assemble the problem */
timecounter tc;
tc.tic();

auto assembler = make_obstacle_assembler(msh, in_A, hdi);



for (auto& cl : msh.cells)
{
auto gr = make_hho_laplacian(msh, cl, hdi);
Matrix<T, Dynamic, Dynamic> stab = make_hho_fancy_stabilization(msh, cl, gr.first, hdi);
Matrix<T, Dynamic, Dynamic> lc = gr.second + stab;
Matrix<T, Dynamic, 1> f = Matrix<T, Dynamic, 1>::Zero(lc.rows());
f = make_rhs(msh, cl, hdi.cell_degree(), rhs_fun);
assembler.assemble_A(msh, cl, lc, f, gamma, bcs_fun);
f = make_rhs(msh, cl, hdi.cell_degree(), rhs_fun, quadrature_degree_increase);
assembler.assemble(msh, cl, lc, f, gamma, bcs_fun);
}

assembler.finalize();
Expand All @@ -147,16 +160,10 @@ void run_hho_obstacle(const Mesh& msh, size_t degree)
std::cout << bold << yellow << "Assembly: " << tc << " seconds" << reset << std::endl;


ss.str("");
ss << "matrix_" << iter << ".txt";

dump_sparse_matrix(assembler.LHS, ss.str());


silo.add_variable("mesh", "difference", diff.data(), diff.size(), zonal_variable_t);
silo.add_variable("mesh", "active", active.data(), active.size(), zonal_variable_t);


/* Solve it */
tc.tic();

SparseLU<SparseMatrix<T>> solver;
Expand All @@ -168,11 +175,11 @@ void run_hho_obstacle(const Mesh& msh, size_t degree)
tc.toc();
std::cout << bold << yellow << "Solver: " << tc << " seconds" << reset << std::endl;

Matrix<T, Dynamic, Dynamic> alpha_prev = alpha;
/* Postprocess, per iteration part */
Matrix<T, Dynamic, 1> alpha_prev = alpha;
assembler.expand_solution(msh, sol, sol_fun, gamma, alpha, beta);


silo.add_variable("mesh", "alpha", alpha.data(), alpha.size(), zonal_variable_t);
silo.add_variable("mesh", "alpha", alpha.data(), num_cells, zonal_variable_t);
silo.add_variable("mesh", "beta", beta.data(), beta.size(), zonal_variable_t);
silo.add_variable("mesh", "expected_solution", expected_solution.data(), expected_solution.size(), zonal_variable_t);

Expand All @@ -184,7 +191,21 @@ void run_hho_obstacle(const Mesh& msh, size_t degree)
iter++;
}

/* Postprocess, final part */
T error = 0.0;
for (auto& cl : msh.cells)
{
Matrix<T, Dynamic, 1> local = take_local_data(msh, cl, hdi, alpha);
Matrix<T, Dynamic, 1> proj = project_function(msh, cl, hdi, sol_fun, quadrature_degree_increase);
auto gr = make_hho_laplacian(msh, cl, hdi);
Matrix<T, Dynamic, Dynamic> stab = make_hho_fancy_stabilization(msh, cl, gr.first, hdi);
Matrix<T, Dynamic, Dynamic> lc = gr.second + stab;
Matrix<T, Dynamic, 1> diff = local - proj;

error += diff.dot(lc*diff);
}

std::cout << green << "Error: " << std::sqrt(error) << reset << std::endl;



Expand All @@ -203,6 +224,8 @@ int main(int argc, char **argv)
mesh_init_params<T> mip;
mip.Nx = 5;
mip.Ny = 5;
mip.min_x = -1;
mip.min_y = -1;

bool displace_nodes = true;

Expand Down
21 changes: 21 additions & 0 deletions apps/obstacle/results/convergence.plot
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,21 @@
set terminal postscript eps enhanced color font 'Helvetica,16'
set output 'convergence.eps'

set logscale x
set logscale y

set xrange [0.3:0.01] reverse

set format "%g"
#set format "%g"

plot 'convergence.txt' using (2/$1):2 w lp lw 2 title "k=0", \
'convergence.txt' using (2/$1):3 w lp lw 2 title "k=1", \
'-' using 1:2 w lp lw 2 title "ref. 1", \
'-' using 1:2 w lp lw 2 title "ref. 1.5"
0.1 0.8
0.02 0.16
e
0.1 0.02
0.02 0.00179
e
5 changes: 5 additions & 0 deletions apps/obstacle/results/convergence.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
8 2.26205 0.197735
16 1.2833 0.0588187
32 0.650286 0.0171607
64 0.326314 0.00529786
128 0.163344 0.00168321
26 changes: 13 additions & 13 deletions src/core/core_bits/utils.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -112,14 +112,14 @@ class hho_degree_info

template<typename Mesh, typename T = typename Mesh::coordinate_type>
Matrix<T, Dynamic, Dynamic>
make_mass_matrix(const Mesh& msh, const typename Mesh::cell_type& cl, size_t degree)
make_mass_matrix(const Mesh& msh, const typename Mesh::cell_type& cl, size_t degree, size_t di = 0)
{
cell_basis<Mesh,T> cb(msh, cl, degree);
auto cbs = cb.size();

Matrix<T, Dynamic, Dynamic> ret = Matrix<T, Dynamic, Dynamic>::Zero(cbs, cbs);

auto qps = integrate(msh, cl, 2*degree);
auto qps = integrate(msh, cl, 2*(degree+di));

for (auto& qp : qps)
{
Expand All @@ -132,14 +132,14 @@ make_mass_matrix(const Mesh& msh, const typename Mesh::cell_type& cl, size_t deg

template<typename Mesh, typename T = typename Mesh::coordinate_type>
Matrix<T, Dynamic, Dynamic>
make_mass_matrix(const Mesh& msh, const typename Mesh::face_type& fc, size_t degree)
make_mass_matrix(const Mesh& msh, const typename Mesh::face_type& fc, size_t degree, size_t di = 0)
{
face_basis<Mesh,T> fb(msh, fc, degree);
auto fbs = fb.size();

Matrix<T, Dynamic, Dynamic> ret = Matrix<T, Dynamic, Dynamic>::Zero(fbs, fbs);

auto qps = integrate(msh, fc, 2*degree);
auto qps = integrate(msh, fc, 2*(degree+di));

for (auto& qp : qps)
{
Expand All @@ -153,7 +153,7 @@ make_mass_matrix(const Mesh& msh, const typename Mesh::face_type& fc, size_t deg
template<typename Mesh, typename Function>
Matrix<typename Mesh::coordinate_type, Dynamic, 1>
make_rhs(const Mesh& msh, const typename Mesh::cell_type& cl,
size_t degree, const Function& f)
size_t degree, const Function& f, size_t di = 0)
{
using T = typename Mesh::coordinate_type;

Expand All @@ -162,7 +162,7 @@ make_rhs(const Mesh& msh, const typename Mesh::cell_type& cl,

Matrix<T, Dynamic, 1> ret = Matrix<T, Dynamic, 1>::Zero(cbs);

auto qps = integrate(msh, cl, 2*degree);
auto qps = integrate(msh, cl, 2*(degree+di));

for (auto& qp : qps)
{
Expand All @@ -176,7 +176,7 @@ make_rhs(const Mesh& msh, const typename Mesh::cell_type& cl,
template<typename Mesh, typename Function>
Matrix<typename Mesh::coordinate_type, Dynamic, 1>
make_rhs(const Mesh& msh, const typename Mesh::face_type& fc,
size_t degree, const Function& f)
size_t degree, const Function& f, size_t di = 0)
{
using T = typename Mesh::coordinate_type;

Expand All @@ -185,7 +185,7 @@ make_rhs(const Mesh& msh, const typename Mesh::face_type& fc,

Matrix<T, Dynamic, 1> ret = Matrix<T, Dynamic, 1>::Zero(fbs);

auto qps = integrate(msh, fc, 2*degree);
auto qps = integrate(msh, fc, 2*(degree+di));

for (auto& qp : qps)
{
Expand All @@ -199,7 +199,7 @@ make_rhs(const Mesh& msh, const typename Mesh::face_type& fc,
template<typename Mesh, typename Function>
Matrix<typename Mesh::coordinate_type, Dynamic, 1>
project_function(const Mesh& msh, const typename Mesh::cell_type& cl,
hho_degree_info hdi, const Function& f)
hho_degree_info hdi, const Function& f, size_t di = 0)
{
using T = typename Mesh::coordinate_type;

Expand All @@ -208,16 +208,16 @@ project_function(const Mesh& msh, const typename Mesh::cell_type& cl,

Matrix<T, Dynamic, 1> ret = Matrix<T, Dynamic, 1>::Zero(cbs+4*fbs);

Matrix<T, Dynamic, Dynamic> cell_mm = make_mass_matrix(msh, cl, hdi.cell_degree());
Matrix<T, Dynamic, 1> cell_rhs = make_rhs(msh, cl, hdi.cell_degree(), f);
Matrix<T, Dynamic, Dynamic> cell_mm = make_mass_matrix(msh, cl, hdi.cell_degree(), di);
Matrix<T, Dynamic, 1> cell_rhs = make_rhs(msh, cl, hdi.cell_degree(), f, di);
ret.block(0, 0, cbs, 1) = cell_mm.llt().solve(cell_rhs);

auto fcs = faces(msh, cl);
for (size_t i = 0; i < 4; i++)
{
auto fc = fcs[i];
Matrix<T, Dynamic, Dynamic> face_mm = make_mass_matrix(msh, fc, hdi.face_degree());
Matrix<T, Dynamic, 1> face_rhs = make_rhs(msh, fc, hdi.face_degree(), f);
Matrix<T, Dynamic, Dynamic> face_mm = make_mass_matrix(msh, fc, hdi.face_degree(), di);
Matrix<T, Dynamic, 1> face_rhs = make_rhs(msh, fc, hdi.face_degree(), f, di);
ret.block(cbs+i*fbs, 0, fbs, 1) = face_mm.llt().solve(face_rhs);
}

Expand Down
Loading

0 comments on commit 7a071b4

Please sign in to comment.