Skip to content
Snippets Groups Projects
Commit 1ab9f794 authored by Matteo Cicuttin's avatar Matteo Cicuttin
Browse files

Restructured tests.

parent b900d3ed
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
CMakeLists.txt
+ 1
− 1
View file @ 1ab9f794
......@@ -233,7 +233,7 @@ include_directories(${CMAKE_CURRENT_SOURCE_DIR})
add_executable(gmsh_gpu_dg main.cpp gmsh_io.cpp reference_element.cpp physical_element.cpp entity.cpp kernels_cpu.cpp)
target_link_libraries(gmsh_gpu_dg ${LINK_LIBS})
add_executable(test test.cpp gmsh_io.cpp reference_element.cpp physical_element.cpp entity.cpp kernels_cpu.cpp)
add_executable(test test.cpp test_basics.cpp test_mass.cpp test_differentiation.cpp gmsh_io.cpp reference_element.cpp physical_element.cpp entity.cpp kernels_cpu.cpp)
target_link_libraries(test ${LINK_LIBS})
add_executable(dump_orient dump_orient.cpp)
......
kernels_cpu.cpp
+ 1
− 0
View file @ 1ab9f794
#include "types.h"
#include "entity_data.h"
/* Field derivatives kernel, case for elements with planar faces */
void compute_field_derivatives(const entity_data_cpu<1>& ed,
const vecxd& f, vecxd& df_dx, vecxd& df_dy, vecxd& df_dz)
{
......
test.cpp
+ 11
− 304
View file @ 1ab9f794
......@@ -4,341 +4,48 @@
#include "test.h"
#include "gmsh_io.h"
#include "silo_output.hpp"
#include "sgr.hpp"
#include "entity_data.h"
#include "kernels_cpu.h"
using namespace sgr;
#define SPHERE_RADIUS 0.2
static void
resize_mesh(double mesh_h)
{
gvp_t vp;
gm::getEntities(vp);
gmm::setSize(vp, mesh_h);
}
static void
make_geometry(int order, double mesh_h)
{
gm::add("difftest");
std::vector<std::pair<int,int>> objects;
objects.push_back(
std::make_pair(3, gmo::addBox(0.0, 0.0, 0.0, 1.0, 1.0, 1.0) )
);
//objects.push_back(
// std::make_pair(3, gmo::addSphere(0.5, 0.5, 0.5, SPHERE_RADIUS) )
// );
//std::vector<std::pair<int, int>> tools;
//gmsh::vectorpair odt;
//std::vector<gmsh::vectorpair> odtm;
//gmo::fragment(objects, tools, odt, odtm);
gmo::synchronize();
resize_mesh(mesh_h);
}
int test_basics(int geometric_order, int approximation_order)
{
model mod(geometric_order, approximation_order);
auto idx = geometric_order - 1;
std::cout << cyanfg << "Testing geometric order " << geometric_order;
std::cout << ", approximation order = " << approximation_order << nofg;
std::cout << std::endl;
double vol_e0 = (4.0/3.0)*M_PI*SPHERE_RADIUS*SPHERE_RADIUS*SPHERE_RADIUS;
double vol_e1 = 1.0 - vol_e0;
/*****************************************/
std::cout << "Volumes by reference quadrature points and determinants" << std::endl;
std::vector<double> t1_vals { 0.07, 0.00145, 0.00021, 1.83e-5 };
COMPARE_VALUES_RELATIVE("Entity 0 vol r+d", mod.entity_at(0).measure(), vol_e0, t1_vals[idx]);
std::vector<double> t2_vals { 0.024, 5.02e-5, 7.003e-06, 6.35e-7 };
COMPARE_VALUES_RELATIVE("Entity 1 vol r+d", mod.entity_at(1).measure(), vol_e1, t2_vals[idx]);
/*****************************************/
std::cout << "Volumes by physical quadrature points" << std::endl;
double vol_e0_pqp = 0.0;
auto& e0 = mod.entity_at(0);
for (const auto& e : e0)
{
const auto& pqps = e.integration_points();
for (const auto& qp : pqps)
vol_e0_pqp += qp.weight();
}
COMPARE_VALUES_RELATIVE("Entity 0 vol pqp", vol_e0_pqp, vol_e0, t1_vals[idx]);
double vol_e1_pqp = 0.0;
auto& e1 = mod.entity_at(1);
for (const auto& e : e1)
{
const auto& pqps = e.integration_points();
for (const auto& qp : pqps)
vol_e1_pqp += qp.weight();
}
COMPARE_VALUES_RELATIVE("Entity 1 vol pqp", vol_e1_pqp, vol_e1, t2_vals[idx]);
/*****************************************/
std::cout << "Volumes by jacobians" << std::endl;
double vol_e0_j = 0.0;
for (const auto& e : e0)
{
const auto& js = e.jacobians();
auto& re = e0.refelem(e);
auto rqps = re.integration_points();
auto dets = e.determinants();
for (size_t i = 0; i < js.size(); i++)
{
double dj = js[i].determinant();
double dg = dets[i];
auto relerr = dj/dg - 1;
assert( std::abs(relerr) < 3e-14 );
vol_e0_j += dj*rqps[i].weight();
}
}
COMPARE_VALUES_RELATIVE("Entity 0 vol j", vol_e0_j, vol_e0, t1_vals[idx]);
double vol_e1_j = 0.0;
for (const auto& e : e1)
{
const auto& js = e.jacobians();
auto& re = e1.refelem(e);
auto rqps = re.integration_points();
auto dets = e.determinants();
for (size_t i = 0; i < js.size(); i++)
{
double dj = js[i].determinant();
double dg = dets[i];
auto relerr = dj/dg - 1;
assert( std::abs(relerr) < 3e-14 );
vol_e1_j += dj*rqps[i].weight();
}
}
COMPARE_VALUES_RELATIVE("Entity 1 vol j", vol_e1_j, vol_e1, t2_vals[idx]);
/*****************************************/
std::cout << "Volumes by integrating f(x) = 1 over the domain" << std::endl;
double vol_e0_mass = 0.0;
auto f = [](const point_3d& pt) -> double { return 1; };
vecxd f_e0 = e0.project(f);
for (size_t iT = 0; iT < e0.num_elements(); iT++)
{
auto& pe = e0.elem(iT);
auto& re = e0.refelem(pe);
auto num_bf = re.num_basis_functions();
matxd mass = e0.mass_matrix(iT);
vecxd f_pe = f_e0.segment(num_bf*iT, num_bf);
vol_e0_mass += f_pe.transpose() * mass * f_pe;
}
COMPARE_VALUES_RELATIVE("Entity 0 vol mass", vol_e0_mass, vol_e0, t1_vals[idx]);
double vol_e1_mass = 0.0;
vecxd f_e1 = e1.project(f);
for (size_t iT = 0; iT < e1.num_elements(); iT++)
{
auto& pe = e1.elem(iT);
auto& re = e1.refelem(pe);
auto num_bf = re.num_basis_functions();
matxd mass = e1.mass_matrix(iT);
vecxd f_pe = f_e1.segment(num_bf*iT, num_bf);
vol_e1_mass += f_pe.transpose() * mass * f_pe;
}
COMPARE_VALUES_RELATIVE("Entity 1 vol mass", vol_e1_mass, vol_e1, t2_vals[idx]);
return 0;
}
#define OVERINTEGRATE
int test_mass_convergence(int geometric_order, int approximation_order)
{
std::vector<double> sizes({ 0.32, 0.16, 0.08, 0.04 });
std::vector<double> errors;
std::cout << cyanfg << "Testing geometric order " << geometric_order;
std::cout << ", approximation order = " << approximation_order << nofg;
std::cout << std::endl;
auto test_f = [](const point_3d& pt) -> double {
return std::sin(M_PI*pt.x())*std::sin(M_PI*pt.y())*std::sin(M_PI*pt.z());
};
for (size_t i = 0; i < sizes.size(); i++)
{
double h = sizes[i];
make_geometry(0,h);
model mod(geometric_order, approximation_order);
auto& e0 = mod.entity_at(0);
e0.sort_by_orientation();
vecxd test_f_e0 = e0.project(test_f);
double proj_error_e0 = 0.0;
#ifdef OVERINTEGRATE
auto [rps, pqps] = integrate(e0, 2*approximation_order+1);
auto num_qp = rps.size();
#endif
for (size_t iT = 0; iT < e0.num_elements(); iT++)
{
auto& pe = e0.elem(iT);
auto& re = e0.refelem(pe);
auto num_bf = re.num_basis_functions();
vecxd pvals = test_f_e0.segment(iT*num_bf, num_bf);
#ifdef OVERINTEGRATE
for (size_t iQp = 0; iQp < num_qp; iQp++)
{
vecxd phi = re.basis_functions( rps[iQp] );
auto pqpofs = num_qp*iT + iQp;
double diff = (pvals.dot(phi) - test_f(pqps[pqpofs].point()));
proj_error_e0 += std::abs(pqps[pqpofs].weight())*diff*diff;
}
#else
auto pqps = pe.integration_points();
for (size_t iQp = 0; iQp < pqps.size(); iQp++)
{
vecxd phi = re.basis_functions( iQp );
double diff = (pvals.dot(phi) - test_f(pqps[iQp].point()));
proj_error_e0 += std::abs(pqps[iQp].weight())*diff*diff;
}
#endif
}
errors.push_back( std::sqrt(proj_error_e0) );
}
for (auto& error : errors)
std::cout << error << " ";
std::cout << std::endl;
double rate = 0.0;
for (size_t i = 1; i < errors.size(); i++)
std::cout << (rate = std::log2(errors[i-1]/errors[i])) << " ";
std::cout << std::endl;
COMPARE_VALUES_ABSOLUTE("Projection", rate, double(approximation_order+1), 0.2);
return 0;
}
int test_differentiation_convergence(int geometric_order, int approximation_order)
{
std::vector<double> sizes({ 0.32, 0.16, 0.08, 0.04 });
std::vector<double> errors;
std::cout << cyanfg << "Testing geometric order " << geometric_order;
std::cout << ", approximation order = " << approximation_order << nofg;
std::cout << std::endl;
auto f = [](const point_3d& pt) -> double {
return std::sin(M_PI*pt.x())*std::sin(M_PI*pt.y())*std::sin(M_PI*pt.z());
};
auto df_dx = [](const point_3d& pt) -> double {
return M_PI*std::cos(M_PI*pt.x())*std::sin(M_PI*pt.y())*std::sin(M_PI*pt.z());
};
auto df_dy = [](const point_3d& pt) -> double {
return M_PI*std::sin(M_PI*pt.x())*std::cos(M_PI*pt.y())*std::sin(M_PI*pt.z());
};
auto df_dz = [](const point_3d& pt) -> double {
return M_PI*std::sin(M_PI*pt.x())*std::sin(M_PI*pt.y())*std::cos(M_PI*pt.z());
};
for (size_t i = 0; i < sizes.size(); i++)
{
double h = sizes[i];
make_geometry(0,h);
model mod(geometric_order, approximation_order);
auto& e0 = mod.entity_at(0);
e0.sort_by_orientation();
vecxd Pf_e0 = e0.project(f);
vecxd df_dx_e0 = vecxd::Zero( Pf_e0.rows() );
vecxd df_dy_e0 = vecxd::Zero( Pf_e0.rows() );
vecxd df_dz_e0 = vecxd::Zero( Pf_e0.rows() );
entity_data_cpu<1> ed;
e0.populate_entity_data(ed);
compute_field_derivatives(ed, Pf_e0, df_dx_e0, df_dy_e0, df_dz_e0);
vecxd Pdf_dx_e0 = e0.project(df_dx);
vecxd Pdf_dy_e0 = e0.project(df_dy);
vecxd Pdf_dz_e0 = e0.project(df_dz);
double err_x = 0.0;
double err_y = 0.0;
double err_z = 0.0;
for (size_t iT = 0; iT < e0.num_elements(); iT++)
{
auto& pe = e0.elem(iT);
auto& re = e0.refelem(pe);
matxd mass = e0.mass_matrix(iT);
auto num_bf = re.num_basis_functions();
vecxd diff_x = Pdf_dx_e0.segment(iT*num_bf, num_bf) - df_dx_e0.segment(iT*num_bf, num_bf);
vecxd diff_y = Pdf_dy_e0.segment(iT*num_bf, num_bf) - df_dy_e0.segment(iT*num_bf, num_bf);
vecxd diff_z = Pdf_dz_e0.segment(iT*num_bf, num_bf) - df_dz_e0.segment(iT*num_bf, num_bf);
err_x += diff_x.dot(mass*diff_x);
err_y += diff_y.dot(mass*diff_y);
err_z += diff_z.dot(mass*diff_z);
}
std::cout << std::sqrt(err_x) << " " << std::sqrt(err_y) << " " << std::sqrt(err_z) << std::endl;
}
using namespace sgr;
return 0;
}
int main(int argc, char **argv)
int main(void)
{
gmsh::initialize();
gmsh::option::setNumber("General.Terminal", 0);
gmsh::option::setNumber("Mesh.Algorithm", 1);
gmsh::option::setNumber("Mesh.Algorithm3D", 1);
int failed_tests = 0;
/*
#if 0
std::cout << Bmagentafg << " *** TESTING: BASIC OPERATIONS ***" << reset << std::endl;
for (size_t go = 1; go < 5; go++)
{
for (size_t ao = go; ao < 5; ao++)
{
make_geometry(0, 0.1);
failed_tests += test_basics(go, ao);
}
}
*/
/*
std::cout << Bmagentafg << " *** TESTING: PROJECTIONS ***" << reset << std::endl;
for (size_t go = 1; go < 5; go++)
{
for (size_t ao = go; ao < 5; ao++)
failed_tests += test_mass_convergence(go, ao);
}
*/
#endif
std::cout << Bmagentafg << " *** TESTING: DIFFERENTIATION ***" << reset << std::endl;
for (size_t ao = 1; ao < 5; ao++)
failed_tests += test_differentiation_convergence(1, ao);
/*****************************************/
//COMPARE_VALUES_RELATIVE("Entity 0 vol mass", vol_e0_mass, vol_e0, t1_vals[order-1]);
#if 0
#endif
gmsh::finalize();
return failed_tests;
}
......
test.h
+ 6
− 0
View file @ 1ab9f794
......@@ -22,3 +22,9 @@
std::cout << ", ref = " << ref << ", tol = " << tol; \
std::cout << ", real_tol = " << real_tol << std::endl; return 1; }
#define SPHERE_RADIUS 0.2
int test_basics(int, int);
int test_mass_convergence(int, int);
int test_differentiation_convergence(int, int);
test_basics.cpp 0 → 100644
+ 151
− 0
View file @ 1ab9f794
#include <iostream>
#include <unistd.h>
#include "test.h"
#include "gmsh_io.h"
#include "sgr.hpp"
using namespace sgr;
static void
make_geometry(int order, double mesh_h)
{
gm::add("test");
std::vector<std::pair<int,int>> objects;
objects.push_back(
std::make_pair(3, gmo::addBox(0.0, 0.0, 0.0, 1.0, 1.0, 1.0) )
);
objects.push_back(
std::make_pair(3, gmo::addSphere(0.5, 0.5, 0.5, SPHERE_RADIUS) )
);
std::vector<std::pair<int, int>> tools;
gmsh::vectorpair odt;
std::vector<gmsh::vectorpair> odtm;
gmo::fragment(objects, tools, odt, odtm);
gmo::synchronize();
gvp_t vp;
gm::getEntities(vp);
gmm::setSize(vp, mesh_h);
}
int
test_basics(int geometric_order, int approximation_order)
{
make_geometry(0, 0.1);
model mod(geometric_order, approximation_order);
auto idx = geometric_order - 1;
std::cout << cyanfg << "Testing geometric order " << geometric_order;
std::cout << ", approximation order = " << approximation_order << nofg;
std::cout << std::endl;
double vol_e0 = (4.0/3.0)*M_PI*SPHERE_RADIUS*SPHERE_RADIUS*SPHERE_RADIUS;
double vol_e1 = 1.0 - vol_e0;
/*****************************************/
std::cout << yellowfg << "Volumes by reference quadrature points and determinants";
std::cout << nofg << std::endl;
std::vector<double> t1_vals { 0.07, 0.00145, 0.00024, 1.83e-5 };
COMPARE_VALUES_RELATIVE("Entity 0 vol r+d", mod.entity_at(0).measure(), vol_e0, t1_vals[idx]);
std::vector<double> t2_vals { 0.0024, 5.02e-5, 8.21e-06, 6.35e-7 };
COMPARE_VALUES_RELATIVE("Entity 1 vol r+d", mod.entity_at(1).measure(), vol_e1, t2_vals[idx]);
/*****************************************/
std::cout << yellowfg << "Volumes by physical quadrature points" << nofg << std::endl;
double vol_e0_pqp = 0.0;
auto& e0 = mod.entity_at(0);
for (const auto& e : e0)
{
const auto& pqps = e.integration_points();
for (const auto& qp : pqps)
vol_e0_pqp += qp.weight();
}
COMPARE_VALUES_RELATIVE("Entity 0 vol pqp", vol_e0_pqp, vol_e0, t1_vals[idx]);
double vol_e1_pqp = 0.0;
auto& e1 = mod.entity_at(1);
for (const auto& e : e1)
{
const auto& pqps = e.integration_points();
for (const auto& qp : pqps)
vol_e1_pqp += qp.weight();
}
COMPARE_VALUES_RELATIVE("Entity 1 vol pqp", vol_e1_pqp, vol_e1, t2_vals[idx]);
/*****************************************/
std::cout << yellowfg << "Volumes by jacobians" << nofg << std::endl;
double vol_e0_j = 0.0;
for (const auto& e : e0)
{
const auto& js = e.jacobians();
auto& re = e0.refelem(e);
auto rqps = re.integration_points();
auto dets = e.determinants();
for (size_t i = 0; i < js.size(); i++)
{
double dj = js[i].determinant();
double dg = dets[i];
auto relerr = dj/dg - 1;
assert( std::abs(relerr) < 3e-14 );
vol_e0_j += dj*rqps[i].weight();
}
}
COMPARE_VALUES_RELATIVE("Entity 0 vol j", vol_e0_j, vol_e0, t1_vals[idx]);
double vol_e1_j = 0.0;
for (const auto& e : e1)
{
const auto& js = e.jacobians();
auto& re = e1.refelem(e);
auto rqps = re.integration_points();
auto dets = e.determinants();
for (size_t i = 0; i < js.size(); i++)
{
double dj = js[i].determinant();
double dg = dets[i];
auto relerr = dj/dg - 1;
assert( std::abs(relerr) < 3e-14 );
vol_e1_j += dj*rqps[i].weight();
}
}
COMPARE_VALUES_RELATIVE("Entity 1 vol j", vol_e1_j, vol_e1, t2_vals[idx]);
/*****************************************/
std::cout << yellowfg << "Volumes by integrating f(x) = 1 over the domain";
std::cout << nofg << std::endl;
double vol_e0_mass = 0.0;
auto f = [](const point_3d& pt) -> double { return 1; };
vecxd f_e0 = e0.project(f);
for (size_t iT = 0; iT < e0.num_elements(); iT++)
{
auto& pe = e0.elem(iT);
auto& re = e0.refelem(pe);
auto num_bf = re.num_basis_functions();
matxd mass = e0.mass_matrix(iT);
vecxd f_pe = f_e0.segment(num_bf*iT, num_bf);
vol_e0_mass += f_pe.transpose() * mass * f_pe;
}
COMPARE_VALUES_RELATIVE("Entity 0 vol mass", vol_e0_mass, vol_e0, t1_vals[idx]);
double vol_e1_mass = 0.0;
vecxd f_e1 = e1.project(f);
for (size_t iT = 0; iT < e1.num_elements(); iT++)
{
auto& pe = e1.elem(iT);
auto& re = e1.refelem(pe);
auto num_bf = re.num_basis_functions();
matxd mass = e1.mass_matrix(iT);
vecxd f_pe = f_e1.segment(num_bf*iT, num_bf);
vol_e1_mass += f_pe.transpose() * mass * f_pe;
}
COMPARE_VALUES_RELATIVE("Entity 1 vol mass", vol_e1_mass, vol_e1, t2_vals[idx]);
return 0;
}
test_differentiation.cpp 0 → 100644
+ 117
− 0
View file @ 1ab9f794
#include <iostream>
#include <iomanip>
#include <unistd.h>
#include "test.h"
#include "gmsh_io.h"
#include "sgr.hpp"
#include "entity_data.h"
#include "kernels_cpu.h"
using namespace sgr;
static void
make_geometry(int order, double mesh_h)
{
gm::add("difftest");
gmo::addBox(0.0, 0.0, 0.0, 1.0, 1.0, 1.0);
gmo::synchronize();
gvp_t vp;
gm::getEntities(vp);
gmm::setSize(vp, mesh_h);
}
int test_differentiation_convergence(int geometric_order, int approximation_order)
{
std::vector<double> sizes({ 0.32, 0.16, 0.08, 0.04 });
std::vector<double> errors;
std::cout << cyanfg << "Testing geometric order " << geometric_order;
std::cout << ", approximation order = " << approximation_order << nofg;
std::cout << std::endl;
auto f = [](const point_3d& pt) -> double {
return std::sin(M_PI*pt.x())*std::sin(M_PI*pt.y())*std::sin(M_PI*pt.z());
};
auto df_dx = [](const point_3d& pt) -> double {
return M_PI*std::cos(M_PI*pt.x())*std::sin(M_PI*pt.y())*std::sin(M_PI*pt.z());
};
auto df_dy = [](const point_3d& pt) -> double {
return M_PI*std::sin(M_PI*pt.x())*std::cos(M_PI*pt.y())*std::sin(M_PI*pt.z());
};
auto df_dz = [](const point_3d& pt) -> double {
return M_PI*std::sin(M_PI*pt.x())*std::sin(M_PI*pt.y())*std::cos(M_PI*pt.z());
};
std::vector<double> errs_x;
std::vector<double> errs_y;
std::vector<double> errs_z;
for (size_t i = 0; i < sizes.size(); i++)
{
double h = sizes[i];
make_geometry(0,h);
model mod(geometric_order, approximation_order);
auto& e0 = mod.entity_at(0);
e0.sort_by_orientation();
vecxd Pf_e0 = e0.project(f);
vecxd df_dx_e0 = vecxd::Zero( Pf_e0.rows() );
vecxd df_dy_e0 = vecxd::Zero( Pf_e0.rows() );
vecxd df_dz_e0 = vecxd::Zero( Pf_e0.rows() );
entity_data_cpu<1> ed;
e0.populate_entity_data(ed);
compute_field_derivatives(ed, Pf_e0, df_dx_e0, df_dy_e0, df_dz_e0);
vecxd Pdf_dx_e0 = e0.project(df_dx);
vecxd Pdf_dy_e0 = e0.project(df_dy);
vecxd Pdf_dz_e0 = e0.project(df_dz);
double err_x = 0.0;
double err_y = 0.0;
double err_z = 0.0;
for (size_t iT = 0; iT < e0.num_elements(); iT++)
{
auto& pe = e0.elem(iT);
auto& re = e0.refelem(pe);
matxd mass = e0.mass_matrix(iT);
auto num_bf = re.num_basis_functions();
vecxd diff_x = Pdf_dx_e0.segment(iT*num_bf, num_bf) - df_dx_e0.segment(iT*num_bf, num_bf);
vecxd diff_y = Pdf_dy_e0.segment(iT*num_bf, num_bf) - df_dy_e0.segment(iT*num_bf, num_bf);
vecxd diff_z = Pdf_dz_e0.segment(iT*num_bf, num_bf) - df_dz_e0.segment(iT*num_bf, num_bf);
err_x += diff_x.dot(mass*diff_x);
err_y += diff_y.dot(mass*diff_y);
err_z += diff_z.dot(mass*diff_z);
}
std::cout << std::sqrt(err_x) << " " << std::sqrt(err_y) << " " << std::sqrt(err_z) << std::endl;
errs_x.push_back( std::sqrt(err_x) );
errs_y.push_back( std::sqrt(err_y) );
errs_z.push_back( std::sqrt(err_z) );
}
double rate_x = 0.0;
double rate_y = 0.0;
double rate_z = 0.0;
std::cout << Byellowfg << "rate df/dx rate df/dy rate df/dz" << reset << std::endl;
for (size_t i = 1; i < sizes.size(); i++)
{
std::cout << (rate_x = std::log2(errs_x[i-1]/errs_x[i]) ) << " ";
std::cout << (rate_y = std::log2(errs_y[i-1]/errs_y[i]) ) << " ";
std::cout << (rate_z = std::log2(errs_z[i-1]/errs_z[i]) ) << std::endl;
}
COMPARE_VALUES_ABSOLUTE("df/dx", rate_x, double(approximation_order), 0.2);
COMPARE_VALUES_ABSOLUTE("df/dy", rate_y, double(approximation_order), 0.2);
COMPARE_VALUES_ABSOLUTE("df/dz", rate_z, double(approximation_order), 0.2);
return 0;
}
\ No newline at end of file
test_mass.cpp 0 → 100644
+ 96
− 0
View file @ 1ab9f794
#include <iostream>
#include <iomanip>
#include <unistd.h>
#include "test.h"
#include "gmsh_io.h"
#include "sgr.hpp"
using namespace sgr;
static void
make_geometry(int order, double mesh_h)
{
gm::add("difftest");
gmo::addBox(0.0, 0.0, 0.0, 1.0, 1.0, 1.0);
gmo::synchronize();
gvp_t vp;
gm::getEntities(vp);
gmm::setSize(vp, mesh_h);
}
#define OVERINTEGRATE
int test_mass_convergence(int geometric_order, int approximation_order)
{
std::vector<double> sizes({ 0.32, 0.16, 0.08, 0.04 });
std::vector<double> errors;
std::cout << cyanfg << "Testing geometric order " << geometric_order;
std::cout << ", approximation order = " << approximation_order << nofg;
std::cout << std::endl;
auto test_f = [](const point_3d& pt) -> double {
return std::sin(M_PI*pt.x())*std::sin(M_PI*pt.y())*std::sin(M_PI*pt.z());
};
for (size_t i = 0; i < sizes.size(); i++)
{
double h = sizes[i];
make_geometry(0,h);
model mod(geometric_order, approximation_order);
auto& e0 = mod.entity_at(0);
e0.sort_by_orientation();
vecxd test_f_e0 = e0.project(test_f);
double proj_error_e0 = 0.0;
#ifdef OVERINTEGRATE
auto [rps, pqps] = integrate(e0, 2*approximation_order+1);
auto num_qp = rps.size();
#endif
for (size_t iT = 0; iT < e0.num_elements(); iT++)
{
auto& pe = e0.elem(iT);
auto& re = e0.refelem(pe);
auto num_bf = re.num_basis_functions();
vecxd pvals = test_f_e0.segment(iT*num_bf, num_bf);
#ifdef OVERINTEGRATE
for (size_t iQp = 0; iQp < num_qp; iQp++)
{
vecxd phi = re.basis_functions( rps[iQp] );
auto pqpofs = num_qp*iT + iQp;
double diff = (pvals.dot(phi) - test_f(pqps[pqpofs].point()));
proj_error_e0 += std::abs(pqps[pqpofs].weight())*diff*diff;
}
#else
auto pqps = pe.integration_points();
for (size_t iQp = 0; iQp < pqps.size(); iQp++)
{
vecxd phi = re.basis_functions( iQp );
double diff = (pvals.dot(phi) - test_f(pqps[iQp].point()));
proj_error_e0 += std::abs(pqps[iQp].weight())*diff*diff;
}
#endif
}
errors.push_back( std::sqrt(proj_error_e0) );
}
std::cout << Byellowfg << "Errors: " << reset;
for (auto& error : errors)
std::cout << error << " ";
std::cout << std::endl;
std::cout << Byellowfg << "Rates: " << reset;
double rate = 0.0;
for (size_t i = 1; i < errors.size(); i++)
std::cout << (rate = std::log2(errors[i-1]/errors[i])) << " ";
std::cout << std::endl;
COMPARE_VALUES_ABSOLUTE("Projection", rate, double(approximation_order+1), 0.2);
return 0;
}
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