phelm.h

Go to the documentation of this file.

#ifndef MKE_H
#define MKE_H
/* -*- charset: utf-8 -*- */
/*$Id: phelm.h,v 2ec1bda0e419 2009/07/13 17:18:25 aozeritsky $*/

#include <stdio.h>

#include "polynom.h"

typedef unsigned int uint;

namespace phelm {

struct Point {
        double x; 
        double y; 

        Point(): x(0), y(0) {}
        Point(double x1, double y1): x(x1), y(y1) {}
        Point(double *x1): x(x1[0]), y(x1[1]) {}

        Point operator / (double k)
        {
                return Point(x / k, y / k);
        }

        Point operator * (double k)
        {
                return Point(x * k, y * k);
        }
};

inline Point operator + (const Point & p1, const Point & p2)
{
        return Point(p1.x + p2.x, p1.y + p2.y);
}

struct MeshPoint {
        std::vector < Point > p;

        MeshPoint() {}
        
        MeshPoint(double x, double y) {
                add(Point(x, y));
        }

        MeshPoint(double *x) {
                add(Point(x));
        }

        void add(const Point & p1) {
                p.push_back(p1);
        }

        double x(int zone = 0) const {
                return p[zone].x;
        }

        double y(int zone = 0) const {
                return p[zone].y;
        }
};

struct Triangle {
        int p[3];  
        int z;     
        std::vector < Polynom > phik; 
        double x[3]; 
        double y[3]; 

        Triangle(int p1, int p2, int p3, int zone = 0)
        {
                p[0] = p1;
                p[1] = p2;
                p[2] = p3;
                z    = zone;
        }

        double X(int i, const std::vector < MeshPoint > & ps) const 
        {
                return ps[p[i]].x(z);
        }

        double Y(int i, const std::vector < MeshPoint > & ps) const 
        {
                return ps[p[i]].y(z);
        }

        void prepare(const std::vector < MeshPoint > & ps)
        {
                std::vector < Polynom > & r = phik;
                r.reserve(3);

                // p0
                r.push_back((P2X - X(1, ps)) * (Y(2, ps) - Y(1, ps)) 
                        - (P2Y - Y(1, ps)) * (X(2, ps) - X(1, ps)));
                // p1
                r.push_back((P2X - X(0, ps)) * (Y(2, ps) - Y(0, ps)) 
                        - (P2Y - Y(0, ps)) * (X(2, ps) - X(0, ps)));
                // p2
                r.push_back((P2X - X(0, ps)) * (Y(1, ps) - Y(0, ps)) 
                        - (P2Y - Y(0, ps)) * (X(1, ps) - X(0, ps)));

                for (uint i = 0; i < 3; ++i)
                {
                        r[i] /= r[i].apply(X(i, ps), Y(i, ps));
                }

                x[0] = X(0, ps); y[0] = Y(0, ps);
                x[1] = X(1, ps); y[1] = Y(1, ps);
                x[2] = X(2, ps); y[2] = Y(2, ps);
        }

        const std::vector < Polynom > & elem1() const
        {
                return phik;
        }

        const Polynom & elem1(int p1) const
        {
                if (p1 == p[0]) {
                        return phik[0];
                } else if (p1 == p[1]) {
                        return phik[1];
                } else {
                        return phik[2];
                }
        }
};

struct Mesh {
        typedef std::vector < Triangle > triangles_t;
        typedef std::vector < MeshPoint > points_t;  
        typedef std::vector < int > points_flags_t;  

        triangles_t tr; 
        points_t ps;    

        points_flags_t ps_flags;

        std::vector < std::vector < int > > adj;
        std::vector < int > inner;
        std::vector < int > outer;
        std::vector < int > p2io;

        bool load(FILE * f);

        void prepare();

        void info();
};
 /* main */

typedef double (* x_t)(double u, double v);

void print_function(FILE * to, double * ans, const Mesh & m, 
                                        x_t x = 0, x_t y = 0, x_t z = 0);

void print_function(const char * fname, double * ans, const Mesh & m, 
                                        x_t x = 0, x_t y = 0, x_t z = 0);

void print_inner_function(FILE * to, double * ans, const Mesh & m, 
                                        x_t x = 0, x_t y = 0, x_t z = 0);

void print_inner_function(const char * to, double * ans, const Mesh & m,
                x_t x = 0, x_t y = 0, x_t z = 0);
 /* print */

double generic_scalar_cb(const Polynom & phi_i, const Polynom & phi_j,
                                                 const Triangle & trk, const Mesh & m, int i1, int j1,
                                                 int i2, int j2, void * user_data);

double sphere_scalar_cb(const Polynom & phi_i, const Polynom & phi_j,
                                                const Triangle & trk, const Mesh & m,
                                                int i1, int j1, int i2, int j2, void * user_data);

class Matrix;

double fast_scalar(const double * u, const double * v,
                                   const Mesh & m, Matrix & mat);

double fast_norm(const double * u, const Mesh & m, Matrix & mat);

double fast_dist(const double * u, const double * v,
                                 const Mesh & m, Matrix & mat);
 /* scalar */

typedef double (* f_xy_t)(double x, double y);

typedef double (* f_xyt_t)(double x, double y, double t);

void p2u(double * u, const double * p, const double * bnd, const Mesh & m);

void u2p(double * p, const double * u, const Mesh & m);

void proj(double * F, const Mesh & mesh, f_xy_t f);

void proj_bnd(double * F, const Mesh & m, f_xy_t f);

void proj_bnd(double * F, const double * F1, const Mesh & m);

void set_bnd(double * F, const double * bnd, const Mesh & m);

void proj(double * F, const Mesh & mesh, f_xyt_t f, double t);

void proj_bnd(double * F, const Mesh & m, f_xyt_t f, double t);
 /* proj */

}

#include "phelm_generators.h"

#endif /* MKE_H */

---

Phelm Library
Copyright © 2009 Alexey Ozeritsky. All rights reserved.
http://resetius.ru