You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
100 lines
3.3 KiB
100 lines
3.3 KiB
9 months ago
|
// Provides a very simple demonstration of how to reinitialise a level set function with
|
||
|
// Algoim's high-order closest point algorithms. The file contains a single main() routine;
|
||
|
// compile it as you would for any .cpp file with a main() entry point.
|
||
|
|
||
|
#include <vector>
|
||
|
#include <algorithm>
|
||
|
#include "hocp.hpp"
|
||
|
|
||
|
// Two-dimensional example using a circle of radius 1
|
||
|
void test2d(int n)
|
||
|
{
|
||
|
std::vector<double> phi_data(n*n), sdfexact_data(n*n);
|
||
|
auto phi = [&](const algoim::uvector<int,2>& i) -> double& { return phi_data[i(0)*n + i(1)]; };
|
||
|
auto sdfexact = [&](const algoim::uvector<int,2>& i) -> double& { return sdfexact_data[i(0)*n + i(1)]; };
|
||
|
|
||
|
double dx = 4.0 / n;
|
||
|
for (int i = 0; i < n; ++i)
|
||
|
{
|
||
|
double x = -2.0 + (i + 0.5) * dx;
|
||
|
for (int j = 0; j < n; ++j)
|
||
|
{
|
||
|
double y = -2.0 + (j + 0.5) * dx;
|
||
|
phi(algoim::uvector<int,2>{i,j}) = exp(x*x + y*y - 1.0) - 1.0;
|
||
|
sdfexact(algoim::uvector<int,2>{i,j}) = sqrt(x*x + y*y) - 1.0;
|
||
|
}
|
||
|
}
|
||
|
|
||
|
// Different polynomial interpolants can be used
|
||
|
//algoim::reinit<2,-1>(phi, n, dx, 10.0);
|
||
|
//algoim::reinit<2,2>(phi, n, dx, 10.0);
|
||
|
algoim::reinit<2,3>(phi, n, dx, 10.0);
|
||
|
//algoim::reinit<2,4>(phi, n, dx, 10.0);
|
||
|
//algoim::reinit<2,5>(phi, n, dx, 10.0);
|
||
|
|
||
|
double maxerr = 0.0;
|
||
|
for (int i = 0; i < n*n; ++i)
|
||
|
maxerr = std::max(maxerr, std::abs(phi_data[i] - sdfexact_data[i]));
|
||
|
|
||
|
std::cout << "Infinity norm of error = " << maxerr << std::endl;
|
||
|
}
|
||
|
|
||
|
// Three-dimensional example using a sphere of radius 1
|
||
|
void test3d(int n)
|
||
|
{
|
||
|
std::vector<double> phi_data(n*n*n), sdfexact_data(n*n*n);
|
||
|
auto phi = [&](const algoim::uvector<int,3>& i) -> double& { return phi_data[i(0)*n*n + i(1)*n + i(2)]; };
|
||
|
auto sdfexact = [&](const algoim::uvector<int,3>& i) -> double& { return sdfexact_data[i(0)*n*n + i(1)*n + i(2)]; };
|
||
|
|
||
|
double dx = 4.0 / n;
|
||
|
for (int i = 0; i < n; ++i)
|
||
|
{
|
||
|
double x = -2.0 + (i + 0.5) * dx;
|
||
|
for (int j = 0; j < n; ++j)
|
||
|
{
|
||
|
double y = -2.0 + (j + 0.5) * dx;
|
||
|
for (int k = 0; k < n; ++k)
|
||
|
{
|
||
|
double z = -2.0 + (k + 0.5) * dx;
|
||
|
phi(algoim::uvector<int,3>{i,j,k}) = exp(x*x + y*y + z*z - 1.0) - 1.0;
|
||
|
sdfexact(algoim::uvector<int,3>{i,j,k}) = sqrt(x*x + y*y + z*z) - 1.0;
|
||
|
}
|
||
|
}
|
||
|
}
|
||
|
|
||
|
// Different polynomial interpolants can be used
|
||
|
//algoim::reinit<3,-1>(phi, n, dx, 10.0);
|
||
|
//algoim::reinit<3,2>(phi, n, dx, 10.0);
|
||
|
algoim::reinit<3,3>(phi, n, dx, 10.0);
|
||
|
//algoim::reinit<3,4>(phi, n, dx, 10.0);
|
||
|
//algoim::reinit<3,5>(phi, n, dx, 10.0);
|
||
|
|
||
|
double maxerr = 0.0;
|
||
|
for (int i = 0; i < n*n*n; ++i)
|
||
|
maxerr = std::max(maxerr, std::abs(phi_data[i] - sdfexact_data[i]));
|
||
|
|
||
|
std::cout << "Infinity norm of error = " << maxerr << std::endl;
|
||
|
}
|
||
|
|
||
|
#if ALGOIM_EXAMPLES_DRIVER == 0 || ALGOIM_EXAMPLES_DRIVER == 2
|
||
|
// A very simple driver
|
||
|
int main(int argc, char* argv[])
|
||
|
{
|
||
|
while (true)
|
||
|
{
|
||
|
int dim = 0, n = 0;
|
||
|
std::cout << "Enter the dimension N = 2 or N = 3: ";
|
||
|
std::cin >> dim;
|
||
|
if (dim != 2 && dim != 3) break;
|
||
|
std::cout << "Enter number of grid points in each dimension: ";
|
||
|
std::cin >> n;
|
||
|
if (n < 2) break;
|
||
|
if (dim == 2)
|
||
|
test2d(n);
|
||
|
else
|
||
|
test3d(n);
|
||
|
}
|
||
|
return 0;
|
||
|
}
|
||
|
#endif
|