zhangyiss/magtet
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

add sample code m3d

Browse Source
This commit is contained in:
张壹
2021-07-28 07:54:12 +08:00
parent d645456208
commit 5b65eab987
36 changed files with 597430 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

166
m3d/src/grid.cpp Executable file
View File

@@ -0,0 +1,166 @@
#include "grid.h"
Grid::Grid (std::string node_file, std::string tet_file):
_node_file (node_file),
_tet_file (tet_file)
{
// do the job
this->read_node();
this->read_tet();
}
Grid::~Grid ()
{
// delete tets
for(unsigned int i=0;i<(_elements.size());++i)
if(_elements[i]!=NULL) {
delete _elements[i]; _elements[i]=NULL;
}
this->_elements.clear();
// delete the nodes.
for(unsigned int i=0;i<(_nodes.size());++i) {
if(_nodes[i]!=NULL) {
delete _nodes[i]; _nodes[i]=NULL;
}
}
this->_nodes.clear();
return;
}
/*
get the i-th element
*/
Tet& Grid::get_elem (const unsigned int i)
{
assert (i < this->n_elems());
assert (_elements[i] != NULL);
return *this->_elements[i];
}
/*
read nodes or sites
*/
void Grid::read_node()
{
// Check input buffer
std::ifstream node_stream(this->_node_file.c_str());
assert(node_stream.good());
unsigned int n_nodes=0;
node_stream >> n_nodes;
assert(n_nodes>0);
this->_nodes.clear();
this->_nodes.resize(n_nodes);
for (unsigned int i=0; i<n_nodes; i++) {
unsigned int node_lab=0;
Point xyz;
node_stream >> node_lab
>> xyz(0) // x-coordinate value
>> xyz(1) // y-coordinate value
>> xyz(2); // z-coordinate value
Node* newnode = new Node(xyz, node_lab);
this->_nodes[i]= newnode;
}
return;
}
/*
read elements and magnetization vector M
*/
void Grid::read_tet ()
{
// Check input buffer
std::ifstream ele_stream(this->_tet_file.c_str());
assert (ele_stream.good());
assert (this->_nodes.size()>0);
unsigned int n_elems=0;
ele_stream >> n_elems;
assert(n_elems>0);
this->_elements.clear();
this->_elements.resize(n_elems);
for(unsigned int i=0; i<n_elems; i++) {
unsigned int element_id;
Node* v[4];
ele_stream >> element_id;
for (unsigned int j=0; j<4; j++) {
unsigned int node_id;
ele_stream >> node_id;
v[j] = this->_nodes[node_id];
}
double mx, my, mz;
ele_stream >> mx >> my >> mz;
Tet* tet = new Tet(element_id, v, Point(mx,my,mz));
this->_elements[i] = tet;
}
this->_number_elems = n_elems;
return;
}
std::ostream& operator <<(std::ostream& os, Grid& mesh)
{
os << "n_nodes()=" << mesh.n_nodes()
<< " \nn_elems()=" << mesh.n_elems()
<< " \n";
return os;
}
/*
Write the tetrahedral elements into VTk format so that
the grid can be shown by Paraview softeware.
*/
void Grid::write_out_vtk(const std::string name)
{
// Open the of stream
std::ofstream vtk_mesh((name+".vtk").c_str());
if(!vtk_mesh.good()) {
std::cerr<<"Can not open file:\t"<<name+".vtk"
<<std::endl;
} else {
vtk_mesh<<"# vtk DataFile Version 3.0\n" //File version and identifier
<<"Zhengyong Ren\n"
<<"ASCII\n";
vtk_mesh<<"DATASET UNSTRUCTURED_GRID\n";
vtk_mesh<<"\nPOINTS\t"<<this->_nodes.size()<<"\tdouble\n";
for(unsigned int i=0; i<this->_nodes.size(); i++) {
const Node& temp_node= (*this->_nodes[i]);
vtk_mesh<<temp_node(0)<<"\t"
<<temp_node(1)<<"\t"
<<temp_node(2)<<"\n";
}
vtk_mesh<<"\nCELLS\t"<<this->n_elems()<<"\t"<<this->n_elems()*5<<"\n";
for(unsigned int i=0; i<this->n_elems(); i++) {
Tet& temp_elem= this->get_elem(i);
vtk_mesh<<(unsigned int)4<<"\t"
<<temp_elem.get_node_id(0)<<"\t"
<<temp_elem.get_node_id(1)<<"\t"
<<temp_elem.get_node_id(2)<<"\t"
<<temp_elem.get_node_id(3)<<"\n";
}
vtk_mesh<<"\nCELL_TYPES\t"<<this->n_elems()<<"\n";
for(unsigned int i=0; i<this->n_elems(); i++) {
vtk_mesh<<(unsigned int)10<<"\n"; // 10--tet
}
vtk_mesh<<"\n";
}
vtk_mesh.close();
// done
}

37
m3d/src/grid.h Executable file
View File

@@ -0,0 +1,37 @@
#ifndef _GRID_H
#define _GRID_H
#include <vector>
#include <string>
#include <cassert>
#include <iostream>
#include <fstream>
#include "node.h"
#include "tet.h"
class Grid
{
public:
Grid(std::string node_file, std::string tet_file);
~Grid();
// tets
unsigned int n_nodes() { return _nodes.size(); }
unsigned int n_elems() { return _number_elems; }
Tet& get_elem (unsigned int i);
void read_node();
void read_tet ();
friend std::ostream& operator <<(std::ostream& os, const Grid& mesh);
void write_out_vtk(const std::string name);
public:
std::string _node_file;
std::string _tet_file;
std::vector<Node*> _nodes;
std::vector<Tet*> _elements;
unsigned int _number_elems;
};
#endif // _GRID_H

365
m3d/src/m3d.cpp Executable file
View File

@@ -0,0 +1,365 @@
// C++ includes
#include "m3d.h"
#include <iostream>
/* m3d
* These files may be found at:
* http://software.seg.org/2017/0006
*
* Please first read and accept:
* http://software.seg.org/disclaimer.txt
*/
/*
Input: a tetrahedon and an observation site
Output: ni (outgoing normal vector of 4 surfaces)
wi0 (heights from site to 4 surfaces
Rij_minus and Rij_plus
oi (projection point of site on 4 surfaces)
Sij_minus and Sij_plus
Rij0 (distance from site to edges of 4 surfaces)
tij (tangential normal vector on edges)
mij (outgoing normal vector of edges)
*/
M3D::M3D(Tet& tet_, Point r_):
tet (tet_),
r (r_),
TOL (1e-10),
Cm (1e-07)
{
// calculate variables for facets and edges
for(int i=0; i<4; i++) {
Point n;
this->tet.get_facet_normal_vector(n, i);
ni[i] = n;
double w0;
for(int j=0; j<3; j++) {
Node* a[2];
this->tet.get_facet_edge_ends (a, i, j);
Rij_minus[i][j] = (this->r-(*a[0])).size();
Rij_plus[i][j] = (this->r-(*a[1])).size();
if(j==0) { w0 = (r-(*a[0]))*n;}
}
wi0[i]=w0;
Point o = r - n*w0;
oi[i] = o;
for(int j=0; j<3; j++) {
Point m, t; Node* a[2];
this->tet.get_facet_edge_normal(m, i, j);
this->tet.get_facet_edge_tangential(t, i, j);
this->tet.get_facet_edge_ends (a, i, j);
Sij_minus[i][j] = ((*a[0])-o)*t;
Sij_plus[i][j] = ((*a[1])-o)*t;
mij0[i][j]= ((*a[0])-o)*m;
Rij0[i][j] = std::sqrt(w0*w0+mij0[i][j]*mij0[i][j]);
mij[i][j] = m;
tij[i][j] = t;
}
}
}
/*
Input: surface index i
Output: k0=\iint_{T_{i}} 1/R which is given in equations (24) and (25)
*/
void M3D::compute_1_R(unsigned int i, double &k0)
{
double w0 = wi0[i];
double temp_k=0.;
for(int j=0; j<3; j++) {
double part1=0, part2=0;
if(std::abs(w0)>TOL) {
double beta = std::atan((mij0[i][j]*Sij_plus[i][j])/(Rij0[i][j]*Rij0[i][j]+std::abs(w0)*Rij_plus[i][j]))
- std::atan((mij0[i][j]*Sij_minus[i][j])/(Rij0[i][j]*Rij0[i][j]+std::abs(w0)*Rij_minus[i][j]));
part2 = std::abs(w0)*beta;
}
double mj0 = mij0[i][j];
if(std::abs(mj0)>TOL) {
part1 = mj0*std::log((Rij_plus[i][j]+Sij_plus[i][j])/(Rij_minus[i][j]+Sij_minus[i][j]));
}
temp_k += part1 - part2;
}
k0 = temp_k;
return;
}
/*
Input: surface index i
Output: k3=\iint_{T_{i}} 1/R^3 which is given in equations (10), (17), (18) and (19)
*/
void M3D::compute_1_R3(unsigned int i, double &k3)
{
double w0 = wi0[i];
double w0_not_zero_k3=0.;
double w0_is_zero_k3 =0.;
// 1_R3!=0
if(std::abs(w0)>TOL) {
for(int j=0; j<3; j++) {
double beta = std::atan((mij0[i][j]*Sij_plus[i][j])/(Rij0[i][j]*Rij0[i][j]+std::abs(w0)*Rij_plus[i][j]))
- std::atan((mij0[i][j]*Sij_minus[i][j])/(Rij0[i][j]*Rij0[i][j]+std::abs(w0)*Rij_minus[i][j]));
w0_not_zero_k3 += (1.0/std::abs(w0))*beta;
}
}
if(std::abs(w0)<TOL) {
// when w0==0, site r must be outside the triangle
// otherwise, there is a singularity in calculating 1/R3 over triangle
for(int j=0; j<3; j++) {
double m0 = mij0[i][j];
if(std::abs(m0)>TOL) {
double temp = (Sij_plus[i][j]/Rij_plus[i][j]-Sij_minus[i][j]/Rij_minus[i][j]);
w0_is_zero_k3 += temp*(-1.0/m0);
}
}
}
if(std::abs(w0)>TOL) k3 = w0_not_zero_k3;
else k3 = w0_is_zero_k3;
return;
}
/*
Input: surface index i
Output: Ki=\iint_{T_{i}} \nabla 1/ R which is given in equations (5)
where \nabla is on the observation site r
*/
void M3D::compute_grad_1_R(unsigned int i, Point &Ki)
{
Point n = ni[i];
double w0 = wi0[i];
Point temp_k_minus;
for(int j=0; j<3; j++) {
Point part1, part2;
if(std::abs(w0)>TOL) {
// betaij see equation (11)
double betaij= atan((mij0[i][j]*Sij_plus[i][j])/(Rij0[i][j]*Rij0[i][j]+std::abs(w0)*Rij_plus[i][j]))
- atan((mij0[i][j]*Sij_minus[i][j])/(Rij0[i][j]*Rij0[i][j]+std::abs(w0)*Rij_minus[i][j]));
part2 = n*sgn(w0)*betaij;
}
double Aij=0.; // see equations (20-22)
if(std::abs(Rij0[i][j])>TOL) {
Aij = std::log((long double)(Rij_plus[i][j]+Sij_plus[i][j])) - std::log((long double)(Rij_minus[i][j]+Sij_minus[i][j]));
}else if(Sij_plus[i][j]>0.&&Sij_minus[i][j]>0.) {
Aij = std::log(Sij_plus[i][j]) - std::log(Sij_minus[i][j]);
}else if(Sij_plus[i][j]<0.&&Sij_minus[i][j]<0.) {
Aij = std::log(Sij_minus[i][j]*-1.0)- std::log(Sij_plus[i][j]*-1.0);
}else {
std::cout<<"site must not on edge!\n";
std::abort();
}
part1 = mij[i][j]*Aij;
temp_k_minus = temp_k_minus + part1 + part2;
}
Ki = temp_k_minus*-1.0;
return;
}
/*
Input: 4 nodes in vector T
Output: volume(size) of tetrahedron
*/
double M3D::get_tri_size(std::vector<Point>& T)
{
assert(T.size()==3);
Point a=T[1]+(T[0]*-1.0);
Point b=T[2]+(T[0]*-1.0);
Point c=a.cross(b);
return c.size()*0.5;
}
/*
Input: surface index i
Output: k2=\iint_{T_{i}} \nabla \nabla 1/ R which is given in equations (25)
where \nabla is on the observation site r
*/
void M3D::compute_grad_grad_1_R(unsigned int i, double k2[3][3])
{
this->set(k2, 0.);
double temp_k2[3][3];
this->set(temp_k2, 0.);
double w0 = wi0[i];
Point n = ni[i];
for(int j=0; j<3; j++) {
Point grad_Aij(0,0,0); // see equations (37) and (38)
if(std::abs(Rij0[i][j])>TOL) {
double factor_n_mij = (Sij_plus[i][j]/(Rij0[i][j]*Rij0[i][j]*Rij_plus[i][j]) -
Sij_minus[i][j]/(Rij0[i][j]*Rij0[i][j]*Rij_minus[i][j]))*-1.0;
double factor_tij = (1.0/Rij_plus[i][j]-1.0/Rij_minus[i][j])*-1.0;
grad_Aij = n*(w0*factor_n_mij) + tij[i][j]*factor_tij - mij[i][j]*(mij0[i][j]*factor_n_mij);
}else {
double factor_tij = (-1.0/Rij_plus[i][j]+1.0/Rij_minus[i][j]);
grad_Aij = tij[i][j]*factor_tij;
}
Point m = mij[i][j];
for(int l=0; l<3; l++) {
for(int k=0; k<3; k++) {
double temp1 = grad_Aij(l)*m(k)*-1.0;
temp_k2[l][k] += temp1;
}
}
if(std::abs(w0)>TOL) {
Node* a[2];
tet.get_facet_edge_ends (a, i, j);
Point grad_Rij_plus = (r - (*a[1]))*(1.0/Rij_plus[i][j]);
Point grad_Rij_minus = (r - (*a[0]))*(1.0/Rij_minus[i][j]);
Point grad_Sij_plus = tij[i][j]*-1.0;
Point grad_Sij_minus = tij[i][j]*-1.0;
Point grad_w0 = this->ni[i];
Point grad_mij0 = this->mij[i][j]*-1.0;
Point grad_Rij0 = (n*w0-mij[i][j]*mij0[i][j])*(1.0/Rij0[i][j]);
Point grad_abs_wi0 = n*this->sgn(w0);
double a_plus = Rij0[i][j]*Rij0[i][j]+std::abs(w0)*Rij_plus[i][j];
double b_plus = mij0[i][j]*Sij_plus[i][j];
Point grad_a_plus = grad_Rij0*(Rij0[i][j]*2.0) + grad_abs_wi0*Rij_plus[i][j]+grad_Rij_plus*std::abs(w0);
Point grad_b_plus = grad_mij0*Sij_plus[i][j]+grad_Sij_plus*mij0[i][j];
double a_minus = Rij0[i][j]*Rij0[i][j]+std::abs(w0)*Rij_minus[i][j];
double b_minus = mij0[i][j]*Sij_minus[i][j];
Point grad_a_minus = grad_Rij0*(Rij0[i][j]*2.0) + grad_abs_wi0*Rij_minus[i][j]+grad_Rij_minus*std::abs(w0);
Point grad_b_minus = grad_mij0*Sij_minus[i][j]+grad_Sij_minus*mij0[i][j];
Point grad_betaij_plus = (grad_b_plus*a_plus-grad_a_plus*b_plus)*(1.0/(a_plus*a_plus+b_plus*b_plus));
Point grad_betaij_minus = (grad_b_minus*a_minus-grad_a_minus*b_minus)*(1.0/(a_minus*a_minus+b_minus*b_minus));
Point grad_betaij = grad_betaij_plus - grad_betaij_minus; // see equation (11)
for(int l=0; l<3; l++)
for(int k=0; k<3; k++) {
double temp2 = this->sgn(w0)*grad_betaij(l)*n(k)*-1.0;
temp_k2[l][k] += temp2;
}
}
} // loop each facets ends
if(std::abs(w0)<TOL) {
double k3 =0.;
this->compute_1_R3(i, k3);
for(int l=0; l<3; l++) {
for(int k=0; k<3; k++) {
double temp3 = k3*n(l)*n(k)*-1.0; // only need for w0=0
temp_k2[l][k]+= temp3;
}
}
}
// copy k2 = temp_k2
for(int l=0; l<3; l++)
for(int k=0; k<3; k++)
k2[l][k] = temp_k2[l][k];
return;
}
/*
a=b
*/
void M3D::set(double a[3][3], double b)
{
for(int i=0; i<3; i++)
for(int j=0; j<3; j++) a[i][j]=b;
return;
}
/*
sign of real number m
*/
double M3D::sgn(const double m)
{
if (std::abs(m)<TOL) return 0.;
else return m/std::abs(m);
}
/*
calculate the magnetic potential V in equation (24)
*/
void M3D::V(double& v)
{
double temp_V=0.;
Point M = this->tet.get_M();
for(int i=0; i<4; i++) {
double k0;
this->compute_1_R(i, k0);
double mdotn = M*ni[i];
temp_V+= k0*mdotn;
}
v = temp_V*Cm;
return;
}
/*
calculate the magnetic field B in equation (4)
*/
void M3D::B(Point& b)
{
Point temp_B(0,0,0);
Point M = this->tet.get_M();
for(int i=0; i<4; i++) {
Point k1(0,0,0);
this->compute_grad_1_R(i, k1);
double mdotn = M*ni[i];
temp_B = temp_B + k1*mdotn;
}
b = temp_B*(Cm*-1.0);
return ;
}
/*
calculate the magnetic gradient tensor T in equation (26)
*/
void M3D::T(double tb[3][3])
{
double temp_tb[3][3];
this->set(temp_tb, 0);
Point M = this->tet.get_M();
for(int i=0; i<4; i++) {
double k2[3][3];
this->compute_grad_grad_1_R(i, k2);
double mdotn = M*ni[i];
for(int l=0; l<3; l++) {
for(int k=0; k<3; k++) {
temp_tb[l][k] += k2[l][k]*mdotn*Cm*-1.0;
}
}
}
// tb = temp_tb
for(int i=0; i<3; i++)
for(int j=0; j<3; j++)
tb[i][j] = temp_tb[i][j];
return ;
}

58
m3d/src/m3d.h Executable file
View File

@@ -0,0 +1,58 @@
#ifndef _M3D_H
#define _M3D_H
// C++ includes
#include <vector>
#include <cassert>
#include <cmath>
#include <iostream>
#include <fstream>
#include <string>
#include <iomanip>
#include "node.h"
#include "tet.h"
#include "grid.h"
class M3D
{
public:
M3D(Tet& tet, Point r);
void V(double& v);
void B(Point& b);
void T(double tb[3][3]);
double sgn(const double m);
void set(double a[3][3], double b);
void compute_1_R(unsigned int i, double &k0);
void compute_grad_1_R(unsigned int i, Point &k1);
void compute_grad_grad_1_R(unsigned int i, double k2[3][3]);
void compute_1_R3(unsigned int i, double &k3);
const double TOL;
const double Cm ;
double get_tri_size(std::vector<Point>& T);
private:
// variables for edges
double Rij_minus[4][3]; // equation (15)
double Rij_plus [4][3]; // equation (14)
double Sij_minus[4][3]; // equation (8)
double Sij_plus [4][3]; // equation (9)
double Rij0 [4][3]; // equation (13)
double mij0 [4][3]; // mij0 is constant along edge eij
Point mij [4][3]; // mij is the outgoing normal vector of edge eij
Point tij [4][3]; // tij is the tangential normal vector of edge eij
// variables for facets
Point ni [4]; // ni is the outgoing normal vector of surface Ti
double wi0 [4]; // equation (8)
Point oi [4]; // oi is the projection point of site r on surface Ti
// observation site and source element
Tet& tet; // a tetrahedron
Point r; // an observation site
};
#endif // _M3D_H

110
m3d/src/node.h Executable file
View File

@@ -0,0 +1,110 @@
#ifndef _NODE_H
#define _NODE_H
// C++ includes
#include <iomanip>
#include <string>
#include "point.h"
#define INVALID_UNIT
class Node : public Point
{
public:
Node (const double x=.0,
const double y=.0,
const double z=.0,
const unsigned int id=0);
Node (const Node& n);
Node (const Point p,
const unsigned int id =0);
~Node ();
friend std::ostream& operator <<(std::ostream& os, const Node& node);
// Overload operators.
Node& operator= (const Node& node);
bool operator == (const Node& node);
bool operator < (const Node& v) const;
void set_id(const unsigned int id) { _id=id; }
unsigned int get_id() { return _id; }
public:
unsigned int _id;
};
//-----------------------------------------------------------------------
inline
Node::Node (const double x, const double y,const double z,
const unsigned int id) :
Point (x,y,z),
_id (id)
{
}
inline
Node::Node (const Node& n):
Point (n),
_id (n._id)
{
}
inline
Node::Node (const Point p,
const unsigned int id ):
Point(p),
_id (id)
{
}
inline
Node::~Node ()
{
}
inline
Node & Node::operator= (const Node& node)
{
(*this)(0) = node(0);
(*this)(1) = node(1);
(*this)(2) = node(2);
_id = node._id;
return *this;
}
inline
bool Node::operator== (const Node& node)
{
if(_id!=node._id)
return false;
else
{
//compare the coordinates.
return Point(*this)==Point(node);
}
}
inline
bool Node::operator < (const Node& rhs) const
{
// Only compare xyz location.
return (Point)(*this)<(Point)(rhs);
}
inline
std::ostream& operator <<(std::ostream& os, const Node& node)
{
os<<node._id
<<"\t"<<node._coords[0]
<<"\t"<<node._coords[1]
<<"\t"<<node._coords[2];
return os;
}
#endif // NODE_H

218
m3d/src/point.h Executable file
View File

@@ -0,0 +1,218 @@
#ifndef _POINT_H
#define _POINT_H
// C++ includes
#include <cmath>
#include <cassert>
#include <iostream>
class Point
{
public:
Point (const double x=0., const double y=0., const double z=0.);
Point (const Point& p);
virtual ~Point ();
friend std::ostream& operator <<(std::ostream& os, const Point& point);
// Overload operators.
Point& operator=(const Point& p);
double operator () (const unsigned int i) const;
double & operator () (const unsigned int i);
Point operator + (const Point& v) const;
Point operator - (const Point& v) const;
Point operator * (const double a) const;
double operator *(const Point& v) const;
bool operator == (const Point& v) const;
bool operator < (const Point& v) const;
void operator = (const double a);
// Math functions
Point cross(const Point& v) const;
Point unit() const;
double size() const;
void zero();
double* get_xyz() { return _coords; }
protected:
double _coords[3];
const double TOL ;
};
//-----------------------------------------------------------------------
inline
Point::Point (const double x,
const double y,
const double z):
TOL (1e-10)
{
_coords[0] = x; _coords[1] = y; _coords[2] = z;
}
inline
Point::Point (const Point &p):
TOL (1e-10)
{
for (unsigned int i=0; i<3; i++)
_coords[i] = p._coords[i];
}
inline
Point& Point::operator=(const Point& p)
{
_coords[0] = p._coords[0];
_coords[1] = p._coords[1];
_coords[2] = p._coords[2];
return (*this);
}
inline
Point::~Point ()
{
//no space is allocated by the new operator.
//so,there is no delete [].
}
inline
double Point::operator ()(const unsigned int i) const
{
assert (i<3);
return _coords[i];
}
inline
double & Point::operator () (const unsigned int i)
{
assert (i<3);
return _coords[i];
}
inline
Point Point::operator + (const Point &p) const
{
return Point(_coords[0] + p._coords[0],
_coords[1] + p._coords[1],
_coords[2] + p._coords[2]);
}
inline
Point Point::operator - (const Point& p) const
{
return Point(_coords[0] - p._coords[0],
_coords[1] - p._coords[1],
_coords[2] - p._coords[2]);
}
inline
Point Point::operator * (const double factor) const
{
return Point(_coords[0]*factor, _coords[1]*factor,_coords[2]*factor);
}
inline
double Point::operator * (const Point&p) const
{
return (_coords[0]*p(0) + _coords[1]*p(1) + _coords[2]*p(2));
}
inline
bool Point::operator == (const Point& rhs) const
{
return ((fabs(_coords[0] - rhs._coords[0]) +
fabs(_coords[1] - rhs._coords[1]) +
fabs(_coords[2] - rhs._coords[2])) < 3*TOL);
}
inline
bool Point::operator < (const Point& rhs) const
{
//First we assume (this)<rhs true
if (*this == rhs) return false;
if ((*this)(0) < rhs(0)) return true; // <
else if ((*this)(0) > rhs(0)) return false;// >
else if (std::abs((*this)(0)-rhs(0))<1e-6){//vx=rhsx
if ((*this)(1) < rhs(1)) return true;
else if ((*this)(1) > rhs(1)) return false;
else if (std::abs((*this)(1)-rhs(1))<1e-6){//vy=rhsy
if ((*this)(2) < rhs(2)) return true;
else if ((*this)(2) > rhs(2)) return false;
}
}
return false;
}
inline
void Point::operator = (const double a)
{
_coords[0] = a; _coords[1] = a; _coords[2] = a;
}
inline
Point Point::cross(const Point& p) const
{
return Point(_coords[1]*p._coords[2] - _coords[2]*p._coords[1],
-_coords[0]*p._coords[2] + _coords[2]*p._coords[0],
_coords[0]*p._coords[1] - _coords[1]*p._coords[0]);
}
inline
Point Point::unit() const
{
const double length = size();
return Point(_coords[0]/length,
_coords[1]/length,
_coords[2]/length);
}
inline
double Point::size() const
{
double value= std::pow(_coords[0],2)+
std::pow(_coords[1],2)+
std::pow(_coords[2],2);
return std::sqrt(value);
}
inline
void Point::zero()
{
_coords[0]=0.;
_coords[1]=0.;
_coords[2]=0.;
}
inline
std::ostream& operator <<(std::ostream& os, const Point& point)
{
os<<point(0)<<"\t"<<point(1)<<"\t"<<point(2);
return os;
}
#endif // _POINT_H

239
m3d/src/tet.h Executable file
View File

@@ -0,0 +1,239 @@
/*
0
TET
/|\
0 / | \
/ | \2
/ 1| \
1 ..4.|....3
\ | /
\ | /
3 \ | /5
\|/
2
*/
#ifndef _TET_H
#define _TET_H
#include <cstdlib>
#include "node.h"
class Tet
{
public:
Tet (unsigned int id, Node* v[4], Point M);
unsigned int get_id () { return _id;}
Point get_M() { return _M;}
unsigned int get_node_id(unsigned int i);
void get_nodes(Node* nodes[4]);
void get_facet(Node* facet[3], unsigned int i);
void get_facet_normal_vector(Point& n,unsigned int i);
void get_facet_edge_nodes(Node* v[2], unsigned int i, unsigned int j);
void get_facet_edge_normal(Point& m, unsigned int i, unsigned int j);
void get_facet_edge_tangential(Point& t, unsigned int i, unsigned int j);
void get_facet_edge_ends (Node* a[2], unsigned int i, unsigned int j);
double get_size();
private:
unsigned int _id;
Node* _v[4];
Point _M;
};
inline
Tet::Tet(unsigned int id,
Node* v[4],
Point M)
{
this->_id = id;
for(unsigned int i=0; i<4; i++) this->_v[i]=v[i];
this->_M = M;
}
inline double Tet::get_size()
{
const Point & a=(*_v[1])+((*_v[0])*-1.0);
const Point & b=(*_v[1])+((*_v[2])*-1.0);
const Point & c=(*_v[1])+((*_v[3])*-1.0);
return std::abs( a*(c.cross(b)) )/6.;
}
inline
unsigned int Tet::get_node_id(unsigned int i)
{
assert(i<4);
return this->_v[i]->get_id();
}
inline
void Tet::get_nodes(Node* nodes[4])
{
for(int i=0; i<4; i++) nodes[i] = this->_v[i];
return;
}
inline
void Tet::get_facet(Node* facet[3], unsigned int i)
{
assert(i<4);
switch(i) {
case 0:
facet[0] = _v[1]; facet[1] = _v[2]; facet[2] = _v[3];
break;
case 1:
facet[0] = _v[0]; facet[1] = _v[2]; facet[2] = _v[3];
break;
case 2:
facet[0] = _v[0]; facet[1] = _v[1]; facet[2] = _v[3];
break;
case 3:
facet[0] = _v[0]; facet[1] = _v[1]; facet[2] = _v[2];
break;
default:
std::abort();
}
return;
}
inline
void Tet::get_facet_normal_vector(Point& n,unsigned int i)
{
Point normal, facet_center;
assert(i<4);
switch(i) {
case 0:
normal = ((*_v[3])-(*_v[2])).cross((*_v[1])-(*_v[2])).unit();
facet_center = ((*_v[1])+(*_v[2])+(*_v[3]))*(1.0/3.0);
break;
case 1:
normal = ((*_v[0])-(*_v[3])).cross((*_v[2])-(*_v[3])).unit();
facet_center = ((*_v[0])+(*_v[2])+(*_v[3]))*(1.0/3.0);
break;
case 2:
normal = ((*_v[3])-(*_v[0])).cross((*_v[1])-(*_v[0])).unit();
facet_center = ((*_v[0])+(*_v[1])+(*_v[3]))*(1.0/3.0);
break;
case 3:
normal = ((*_v[2])-(*_v[1])).cross((*_v[0])-(*_v[1])).unit();
facet_center = ((*_v[0])+(*_v[1])+(*_v[2]))*(1.0/3.0);
break;
default:
std::abort();
}
Point center = ((*_v[0])+(*_v[1])+(*_v[2])+(*_v[3]))*(0.25);
// we need outgoing normal vector
if(normal*(facet_center-center)<0.0) normal= normal*-1.0;
n = normal;
return;
}
/* 1
0---------2
\ /
\ /
2 \ / 0
\/
1
*/
inline
void Tet::get_facet_edge_nodes(Node* v[2], unsigned int i, unsigned int j)
{
assert(i<4);
assert(j<3);
Node* facet[3];
this->get_facet(facet, i);
switch(j) {
case 0:
v[0] = facet[1]; v[1] = facet[2];
break;
case 1:
v[0] = facet[2]; v[1] = facet[0];
break;
case 2:
v[0] = facet[0]; v[1] = facet[1];
break;
default:
std::abort();
}
return;
}
inline
void Tet::get_facet_edge_normal(Point& m, unsigned int i, unsigned int j)
{
assert(i<4);
assert(j<3);
Node* v[2];
this->get_facet_edge_nodes(v, i, j);
Point t = ((*v[1])-(*v[0])).unit();
Point n;
this->get_facet_normal_vector(n, i);
Point temp_m = t.cross(n);
// we need outgoing normal on edge
Node* facet[3];
this->get_facet(facet, i);
Point facet_center = ((*facet[0])+(*facet[1])+(*facet[2]))*(1.0/3.0);
Point edge_center = ((*v[1])+(*v[0]))*0.5;
if(temp_m*(edge_center-facet_center)<0.0) temp_m= temp_m*-1.0;
m = temp_m;
return;
}
inline
void Tet::get_facet_edge_tangential(Point& t, unsigned int i, unsigned int j)
{
assert(i<4);
assert(j<3);
Point m;
this->get_facet_edge_normal(m, i, j);
Point n;
this->get_facet_normal_vector(n, i);
t = n.cross(m);
}
/*
a[0]-----------a[1]
(a-) (a+)
*/
inline
void Tet::get_facet_edge_ends (Node* a[2], unsigned int i, unsigned int j)
{
assert(i<4);
assert(j<3);
Point t;
this->get_facet_edge_tangential(t, i, j);
Node* v[2];
this->get_facet_edge_nodes(v, i, j);
Point temp_t = ((*v[1])-(*v[0])).unit();
if(t*temp_t>0.0) {
a[0] = v[0];
a[1] = v[1];
}else if(t*temp_t<0.) {
a[0] = v[1];
a[1] = v[0];
}
return;
}
#endif // _Tet_h