test_slaplace.cpp

Laplace equation on a sphere
$\begin{eqnarray*} \Delta \psi &=& f(\varphi, \lambda) \\ \Delta \psi &=& \frac{1}{cos\varphi}\frac{\partial}{\partial\varphi}cos(\varphi)\frac{\partial}{\partial\varphi}\psi+ \frac{1}{cos^2\varphi}\frac{\partial^2}{\partial\lambda^2}\psi\\ \psi|_{\partial\Omega}&=&\psi_0 \\ \end{eqnarray*}$
/* -*- charset: utf-8 -*- */
 
#define _USE_MATH_DEFINES
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include <vector>

#include "phelm.h"
#include "util.h"
#include "laplace.h"

using namespace std;
using namespace phelm;

void usage(const char * name)
{
        fprintf(stderr, "usage: %s [mesh.txt|-]\n", name);
        exit(1);
}

double rp(double x, double y)
{
        //return -6.0 * sin(y) * sin(2.0 * x);
        
        double sx = sin(x);
        double cx = cos(x);
        double sy = sin(y);
        double cy = cos(y);

        double ans = 0.0;

        ans += (sx * sx * (cx * cy - 1) * cy);
        ans += - (cx * cy - 1) * cy;
        ans += -2 * sx * (cx * cy - sx) * (-sx * cy - cx);
        ans += -2 * (cx * cy - sx) * cy;
        if (fabs(ans) > 1e-15) {
                ans /= cx;
        }

        ans += sx * sx * cy * cy;
        ans += -cx * (cx * cy - 1) * cy;
        ans += sy * sy;
        ans += 2 * ipow(-sx * cy - cx, 2);
        ans += 2 * (cx * cy - sx) * (-cx * cy + sx);
        ans += 2 * sy * sy;
        ans += -8 * (sx - 3.0) * sx + 4.0 * cx * cx;
        return ans;
}

double ans(double x, double y)
{
        //return sin(y) * sin(2.0 * x);

        double sx = sin(x);
        double cx = cos(x);
        //double sy = sin(y);
        double cy = cos(y);

        return 0.5 * ipow(cx * cy - 1.0, 2) +
                ipow(cx * cy - sx, 2) +
                2.0 * ipow(sx - 3, 2);
}

double bnd(double x, double y)
{
        return ans(x, y);
}

double nr2(double * a, double * b, int n)
{
        double sum = 0.0;
        for (int i = 0; i < n; ++i) {
                sum += (a[i] - b[i]) * (a[i] - b[i]);
        }
        return sqrt(sum);
}

static double x(double u, double v)
{
        return cos(u) * cos(v);
}

static double y(double u, double v)
{
        return cos(u) * sin(v);
}

static double z(double u, double v)
{
        return sin(u);
}

void test_invert(const Mesh & mesh)
{
        int sz = (int)mesh.ps.size();
        int rs = (int)mesh.inner.size();
        int os = (int)mesh.outer.size();

        vector < double > F(sz);
        vector < double > B(std::max(os, 1));
        vector < double > Ans(sz);
        vector < double > rans(sz);

        proj(&F[0], mesh, rp);
        proj(&rans[0], mesh, ans);
        proj_bnd(&B[0], mesh, bnd);

        Timer t;
        SphereLaplace l(mesh);
        fprintf(stderr, "l->%lf\n", t.elapsed());
        l.solve(&Ans[0], &F[0], &B[0]);

        fprintf(stdout, "L1: invert  err=%.2le\n", dist(&Ans[0], &rans[0], mesh, sphere_scalar_cb, (void*)0));

        {
                FILE * f = fopen("invert_answer.txt", "wb");
                print_function(f, &Ans[0], mesh, x, y, z);
                fclose(f);

                fprintf(stderr, "invert answer saved to 'invert_answer.txt'\n");

                f = fopen("answer.txt", "wb");
                print_function(f, &rans[0], mesh, x, y, z);
                fclose(f);

                f = fopen("rp.txt", "wb");
                print_function(f, &F[0], mesh, x, y, z);
                fclose(f);
        }
}

void test_laplace(const Mesh & mesh)
{
        int sz = (int)mesh.ps.size();
        int rs = (int)mesh.inner.size();
        int os = (int)mesh.outer.size();

        vector < double > U(sz);
        vector < double > LU(sz);
        vector < double > LU1(sz);
        vector < double > B1(std::max(os, 1));
        vector < double > B2(std::max(os, 1));
        vector < double > P(rs);
        vector < double > P1(rs);

        proj(&U[0], mesh, ans);
        proj(&LU[0], mesh, rp);
        proj_bnd(&B1[0], mesh, rp);
        proj_bnd(&B2[0], mesh, ans);

        SphereLaplace l(mesh);
        l.calc1(&LU1[0], &U[0], &B1[0]);

        fprintf(stdout, "L2: laplace err=%.2le\n", dist(&LU[0], &LU1[0], mesh));
        {
                FILE * f = fopen("slu_real.txt", "w");
                print_function(f, &LU[0], mesh, x, y, z);
                fclose(f);
                f = fopen("slu_calc.txt", "w");
                print_function(f, &LU1[0], mesh, x, y, z);
                fclose(f);
        }

        l.solve(&LU[0], &LU1[0], &B2[0]);

        fprintf(stdout, "L3: laplace err=%.2le\n", dist(&U[0], &LU[0], mesh));
}

int main(int argc, char *argv[])
{
        Mesh mesh;
        if (argc > 1) {
                FILE * f = (strcmp(argv[1], "-") == 0) ? stdin : fopen(argv[1], "rb");
                if (!f) {
                        usage(argv[0]);
                }
                if (!mesh.load(f)) {
                        usage(argv[0]);
                }
                fclose(f);
        } else {
                usage(argv[0]);
        }

        test_invert(mesh);
        test_laplace(mesh);

        return 0;
}