From f9bf3b963f07c912088a3433a99f87722db92556 Mon Sep 17 00:00:00 2001
From: Matteo Cicuttin <datafl4sh@toxicnet.eu>
Date: Fri, 23 Apr 2021 22:08:33 +0200
Subject: [PATCH] saving stuff
---
entity.cpp | 31 +++++++++-
physical_element.cpp | 13 ++--
physical_element.h | 1 +
silo_output.hpp | 5 ++
test.cpp | 2 +-
test_lifting.cpp | 144 ++++++++++++++++++++++++++++++++++++++-----
6 files changed, 176 insertions(+), 20 deletions(-)
diff --git a/entity.cpp b/entity.cpp
index 1515448..361ebff 100644
--- a/entity.cpp
+++ b/entity.cpp
@@ -92,6 +92,13 @@ entity::project(const scalar_function& function) const
return ret;
}
+const reference_element&
+entity::cell_refelem(size_t iO) const
+{
+ assert(iO < reference_cells.size());
+ return reference_cells[iO];
+}
+
const reference_element&
entity::cell_refelem(const physical_element& pe) const
{
@@ -458,7 +465,7 @@ entity::populate_lifting_matrices(entity_data_cpu& ed) const
for (size_t i = 0; i < fnt.size(); i++)
for (size_t j = 0; j < cnt.size(); j++)
- if (cnt[j] == fnt[i])
+ if (fnt[i] == cnt[j])
{
assert(f2c.size() == iF*num_face_bf + i);
f2c.push_back(j);
@@ -587,6 +594,28 @@ entity::populate_entity_data(entity_data_cpu& ed) const
mass.inverse();
}
+ ed.fluxdofs_mine.resize(4*ed.num_fluxes*num_cells());
+ for (size_t iT = 0; iT < num_cells(); iT++)
+ {
+ auto& pe = physical_cells[iT];
+ auto cell_keys = pe.bf_keys();
+
+ for (size_t iF = 0; iF < 4; iF++)
+ {
+ auto face_num = 4*iT+iF;
+ auto& pf = physical_faces[face_num];
+ auto face_keys = pf.bf_keys();
+
+ auto face_base = (4*iT + iF) * ed.num_fluxes;
+ auto cell_base = iT * ed.num_bf;
+
+ for (size_t i = 0; i < face_keys.size(); i++)
+ for (size_t j = 0; j < cell_keys.size(); j++)
+ if ( face_keys.at(i) == cell_keys.at(j) )
+ ed.fluxdofs_mine.at(face_base + i) = cell_base + j;
+ }
+ }
+
}
diff --git a/physical_element.cpp b/physical_element.cpp
index 651a94b..5859b58 100644
--- a/physical_element.cpp
+++ b/physical_element.cpp
@@ -108,10 +108,11 @@ physical_element::jacobians(void) const
return m_jacobians;
}
-
-
-
-
+point_3d
+physical_element::barycenter(void) const
+{
+ return m_barycenter;
+}
@@ -127,6 +128,7 @@ physical_elements_factory::physical_elements_factory(const entity_params& ep)
gmm::getBasisFunctionsOrientationForElements(elemType, basis_func_name(approx_order),
cellOrientations, tag);
assert(cellOrientations.size() == cellTags.size());
+ gmm::getBarycenters(elemType, tag, false, false, barycenters);
std::vector<double> coord;
gmm::getKeysForElements(elemType, basis_func_name(approx_order),
@@ -161,6 +163,9 @@ physical_elements_factory::get_elements()
new_pe.m_element_tag = cellTags[elem];
new_pe.m_measure = 0.0;
new_pe.m_gmsh_elemtype = elemType;
+ new_pe.m_barycenter = point_3d(barycenters[3*elem+0],
+ barycenters[3*elem+1],
+ barycenters[3*elem+2]);
new_pe.m_bf_keys.resize(keys_per_elem);
for (size_t i = 0; i < keys_per_elem; i++)
diff --git a/physical_element.h b/physical_element.h
index 1ee451e..4add250 100644
--- a/physical_element.h
+++ b/physical_element.h
@@ -92,6 +92,7 @@ class physical_elements_factory
std::vector<double> cellJacs;
std::vector<double> cellDets;
std::vector<int> cellOrientations;
+ std::vector<double> barycenters;
gmsh::vectorpair keypairs;
int keys_per_elem;
diff --git a/silo_output.hpp b/silo_output.hpp
index daac379..cfa006e 100644
--- a/silo_output.hpp
+++ b/silo_output.hpp
@@ -179,4 +179,9 @@ public:
m_siloDb = nullptr;
}
+
+ ~silo()
+ {
+ close_db();
+ }
};
diff --git a/test.cpp b/test.cpp
index ed5506a..0fc9015 100644
--- a/test.cpp
+++ b/test.cpp
@@ -56,7 +56,7 @@ int main(void)
#endif
std::cout << Bmagentafg << " *** TESTING: LIFTING ***" << reset << std::endl;
- for (size_t ao = 1; ao < 2; ao++)
+ for (size_t ao = 1; ao < 3; ao++)
test_lifting(1, ao);
gmsh::finalize();
diff --git a/test_lifting.cpp b/test_lifting.cpp
index 61158f1..154ccce 100644
--- a/test_lifting.cpp
+++ b/test_lifting.cpp
@@ -1,6 +1,6 @@
#include <iostream>
#include <iomanip>
-
+#include <sstream>
#include <unistd.h>
#include "test.h"
@@ -9,6 +9,7 @@
#include "entity_data.h"
#include "kernels_cpu.h"
#include "timecounter.h"
+#include "silo_output.hpp"
using namespace sgr;
@@ -25,6 +26,39 @@ make_geometry(int order, double mesh_h)
}
+int test_lifting_xx(int geometric_order, int approximation_order)
+{
+ make_geometry(0, 0.25);
+ model mod(geometric_order, approximation_order);
+
+ auto& e0 = mod.entity_at(0);
+
+ entity_data_cpu ed;
+ e0.populate_entity_data(ed);
+
+ for (size_t iO = 0; iO < e0.num_cell_orientations(); iO++)
+ {
+ auto& re = e0.cell_refelem(iO);
+ auto num_bf = re.num_basis_functions();
+
+ auto lift_cols = ed.lifting_matrices.cols();
+
+ matxd mass = re.mass_matrix();
+ matxd curlift = mass * ed.lifting_matrices.block(iO*num_bf, 0, num_bf, lift_cols);
+
+ vecxd ones_cell = vecxd::Ones(num_bf);
+ vecxd ones_face = vecxd::Ones(lift_cols/4);
+
+ for (size_t iF = 0; iF < 4; iF++)
+ {
+ matxd liftblock = curlift.block(0,iF*lift_cols/4,num_bf,lift_cols/4);
+ std::cout << ones_cell.dot(liftblock*ones_face) << std::endl;
+ }
+ }
+
+ return 0;
+}
+
int test_lifting(int geometric_order, int approximation_order)
{
std::vector<double> sizes({ 0.32, 0.16, 0.08, 0.04 });
@@ -34,14 +68,28 @@ int test_lifting(int geometric_order, int approximation_order)
std::cout << ", approximation order = " << approximation_order << nofg;
std::cout << std::endl;
- auto F = [](const point_3d& pt) -> vecxd {
+ /* Test field */
+ auto Fx = [](const point_3d& pt) -> double {
+ return 1; //std::sin(M_PI*pt.x());
+ };
+
+ auto Fy = [](const point_3d& pt) -> double {
+ return 1; //std::sin(M_PI*pt.y());
+ };
+
+ auto Fz = [](const point_3d& pt) -> double {
+ return 1; //std::sin(M_PI*pt.z());
+ };
+
+ auto F = [&](const point_3d& pt) -> vecxd {
vec3d Fret;
- Fret(0) = std::sin(M_PI*pt.y());
- Fret(1) = std::sin(M_PI*pt.z());
- Fret(2) = std::sin(M_PI*pt.x());
+ Fret(0) = Fx(pt);
+ Fret(1) = Fy(pt);
+ Fret(2) = Fz(pt);
return Fret;
};
+ /* Expected divergence */
auto div_F = [](const point_3d& pt) -> double {
return 0;
auto dFx_dx = M_PI*std::cos(M_PI*pt.x());
@@ -50,22 +98,61 @@ int test_lifting(int geometric_order, int approximation_order)
return dFx_dx + dFy_dy + dFz_dz;
};
- for (size_t i = 0; i < sizes.size(); i++)
+ for (size_t sz = 0; sz < sizes.size(); sz++)
{
- double h = sizes[i];
+ double h = sizes[sz];
make_geometry(0,h);
model mod(geometric_order, approximation_order);
+ std::stringstream ss;
+ ss << "lifting_test_" << sz << ".silo";
+
+ silo s;
+ s.create_db( ss.str() );
+ s.import_mesh_from_gmsh();
+ s.write_mesh();
+
auto& e0 = mod.entity_at(0);
- e0.sort_by_orientation();
entity_data_cpu ed;
e0.populate_entity_data(ed);
vecxd exp_field = e0.project(div_F);
- vecxd lift_field = vecxd::Zero( e0.num_cells() * ed.num_bf );
- vecxd flux = vecxd::Zero( e0.num_cells() * 4*ed.num_fluxes );
+ vecxd Fx_proj = e0.project(Fx);
+ vecxd Fy_proj = e0.project(Fy);
+ vecxd Fz_proj = e0.project(Fz);
+
+ auto flx_size = e0.num_cells() * 4*ed.num_fluxes;
+ vecxd Fx_flux = vecxd::Zero( flx_size );
+ vecxd Fy_flux = vecxd::Zero( flx_size );
+ vecxd Fz_flux = vecxd::Zero( flx_size );
+
+ for (size_t i = 0; i < flx_size; i++)
+ {
+ Fx_flux(i) = Fx_proj( ed.fluxdofs_mine[i] );
+ Fy_flux(i) = Fy_proj( ed.fluxdofs_mine[i] );
+ Fz_flux(i) = Fz_proj( ed.fluxdofs_mine[i] );
+ }
+
+ vecxd flux = vecxd::Zero( flx_size );
+ for (size_t iT = 0; iT < e0.num_cells(); iT++)
+ {
+ for (size_t iF = 0; iF < 4; iF++)
+ {
+ vec3d n = ed.normals.row(4*iT+iF);
+ for (size_t k = 0; k < ed.num_fluxes; k++)
+ {
+ auto dof_ofs = (4*iT+iF)*ed.num_fluxes + k;
+ flux(dof_ofs) = Fx_flux(dof_ofs)*n(0) +
+ Fy_flux(dof_ofs)*n(1) +
+ Fz_flux(dof_ofs)*n(2);
+ }
+ }
+ }
+
+ vecxd lift_field = vecxd::Zero( e0.num_cells() * ed.num_bf );
+#if 0
for (size_t iT = 0; iT < e0.num_cells(); iT++)
{
for (size_t iF = 0; iF < 4; iF++)
@@ -94,12 +181,19 @@ int test_lifting(int geometric_order, int approximation_order)
flux.segment(fluxofs, num_fluxes) = mass.ldlt().solve(rhs);
}
}
-
+#endif
timecounter tc;
tc.tic();
compute_flux_lifting(ed, flux, lift_field);
double time = tc.toc();
+ std::vector<double> var_div_F( e0.num_cells() );
+ std::vector<double> var_lift( e0.num_cells() );
+ std::vector<double> var_volp( e0.num_cells() );
+ std::vector<double> var_elift( e0.num_cells() );
+ std::vector<double> var_err( e0.num_cells() );
+ std::vector<double> var_orient( e0.num_cells() );
+
double err = 0.0;
for (size_t iT = 0; iT < e0.num_cells(); iT++)
{
@@ -119,11 +213,33 @@ int test_lifting(int geometric_order, int approximation_order)
vol += pqp.weight() * dphi * F(pqp.point());
}
- vecxd comp_field = lift_field.segment(iT*num_bf, num_bf) - mass.ldlt().solve(vol);
- vecxd diff = comp_field - exp_field.segment(iT*num_bf, num_bf);
- err += diff.dot(mass*diff);
+ vecxd expected = exp_field.segment(iT*num_bf, num_bf);
+ vecxd lf = lift_field.segment(iT*num_bf, num_bf);
+ vecxd volp = mass.ldlt().solve(vol);
+
+ vecxd comp_field = lf - volp;
+ vecxd diff = comp_field - expected;
+ double loc_err = diff.dot(mass*diff);
+ err += loc_err;
+
+ auto bar = pe.barycenter();
+ vecxd phi_bar = re.basis_functions({1./3., 1./3., 1./3.});
+
+ var_div_F.at( pe.original_position() ) = expected.dot(phi_bar);
+ var_lift.at( pe.original_position() ) = lf.dot(phi_bar);
+ var_volp.at( pe.original_position() ) = volp.dot(phi_bar);
+ var_elift.at( pe.original_position() ) = (volp + expected).dot(phi_bar);
+ var_err.at( pe.original_position() ) = loc_err;
+ var_orient.at( pe.original_position() ) = pe.orientation();
}
+ s.write_zonal_variable("div_F", var_div_F);
+ s.write_zonal_variable("lift", var_lift);
+ s.write_zonal_variable("vol", var_volp);
+ s.write_zonal_variable("exp_lift", var_elift);
+ s.write_zonal_variable("err", var_err);
+ s.write_zonal_variable("orient", var_orient);
+
std::cout << std::sqrt(err) << std::endl;
if (geometric_order == 1)
--
GitLab