Commit 55f13dea authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Remove development files

parent 9b623c73
project(Discretization)
set(SOURCE discr_common.cpp
discr_tpfa.cpp
discr_tpfa2.cpp
discr_ntpfa_a.cpp
discr_ntpfa_b.cpp
discr_ntpfa_c.cpp
discr_mpfa_a.cpp
discr_mpfa_b.cpp
discr_nmpfa_a.cpp
discr_nmpfa_b.cpp
discr_mpfa_o.cpp
discr_mpfa_l.cpp
discr_mpfa_c.cpp
discr_triplet_basic.cpp
discr_basic_a.cpp
discr.h)
add_library(Discretization ${SOURCE})
if(USE_MPI)
target_link_libraries(Discretization ${MPI_LIBRARIES})
if(MPI_LINK_FLAGS)
set_target_properties(Discretization PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
\ No newline at end of file
#pragma once
#ifndef _DISCR_H
#include "../../inmost.h"
using namespace INMOST;
//TODO:
//move orphan function inside class
//introduce discr_triplet_basic class
void tensor_prod1(Storage::real * K, Storage::real v[3], Storage::real out[3]);
void tensor_prod2(Storage::real * K, Storage::real v[3], Storage::real out[3]);
void tensor_prod3(Storage::real * K, Storage::real v[3], Storage::real out[3]);
int tensor_prod(int Ktype, Storage::real_array & K,
Storage::real v[3], Storage::real out[3]);
void tensor_add1(Storage::real * K, Storage::real Kout[9]);
void tensor_add2(Storage::real * K, Storage::real Kout[9]);
void tensor_add3(Storage::real * K, Storage::real Kout[9]);
int tensor_add(int Ktype, Storage::real_array & K, Storage::real Kout[9]);
int average_tensor(int Ktype, Storage::real_array & K1,
Storage::real_array & K2, Storage::real Kout[9]);
void vec_diff(Storage::real x[3], Storage::real y[3], Storage::real out[3]);
Storage::real dot_prod(Storage::real x[3], Storage::real y[3]);
//type of the averaging point
#define AVG_NONE 0//material is homogeneous across cells - no point required
#define AVG_CENTROID 1//use centroid
#define AVG_HARMONIC 2//use harmonic averaging point
#define AVG_NONLINEAR 4//use nonlinear flux continuety to define value at the point
#define AVG_REVDIST 8//use reverse distance
#define AVG_NEUMANN 16 //use zero neumann boundary condition
#define CRITICAL_DIST_RATIO 10.0
class discr_basic
{
protected:
Automatizator * aut;
Mesh * m;
Tag K;
Tag bnd_conds;
//Tag AvgCoord, AvgCoefs, AvgStats;
int Ktype;
int error;
ElementType bnd_types;
MIDType bnd_markers, add_markers;
int rem_marker;
int have_bnd;
int allowed_types;
// returns
// 0 if point inside
// 1 if distance to point by distance to farthest node is less then CRITICAL_DIST_RATIO
// 2 otherwise
int isPointInside(Face * f, Storage::real x[3]);
//void InitializeSupportPoints(int allowed_types = AVG_NONE | AVG_CENTROID | AVG_HARMONIC | AVG_NONLINEAR | AVG_REVDIST);
int FindHarmonicPoint(Face * fKL, Cell * cK, Cell * cL,
Tag tensor, int tensor_type,
Storage::real xK[3], Storage::real xL[3],
Storage::real y[3], Storage::real & coef);
int FindBoundaryPoint(Face * fK, Cell * cK,
Tag tensor, int tensor_type,
Storage::real xK[3], Storage::real y[3]);
virtual void Build(Face * f) = 0;
public:
enum DiscrType {TPFA, MPFA, NTPFA, NMPFA};
enum SchemeType{ Linear, Nonlinear };
public:
discr_basic(Automatizator * aut, Mesh * m, std::string tensor_name);
virtual bool needBuild(Face * f) const;
virtual bool needUpdate(Face * f) const;
virtual void Init() = 0;
virtual DiscrType GetDiscrType() const = 0;
virtual SchemeType GetSchemeType() const = 0;
virtual expr Grad(const expr & param) const = 0;
virtual expr BoundaryGrad(const expr & param) const = 0;
virtual expr Interp(ElementType etype,const expr & param) const = 0;
virtual void Update() = 0;
virtual void Export(std::ostream & out, Storage::real trans_scale, Storage::real vol_scale) const; //export to ADGPRS
// by default no boundary information is provided, schemes avoid using boundary
// first marker is the marker that is used inside Automatizator to define dynamic tag domain
// if some boundary element is not used inside scheme it will not be marked by usage_marker
// second marker defines boundary faces/edges, if not set, it is computed automatically
// third parameter defines what kinds of elements may be used
void Boundary(MIDType bnd_marker = 0, ElementType bnd_types = FACE | EDGE);
// mark additional collocation points that should be introduced into the system
MIDType UnknownMarker() const { return add_markers; }
MIDType BoundaryMarker() const { return bnd_markers; }
virtual ~discr_basic();
void AllowTypes(int types) {allowed_types = types;}
};
class discr_triplet_basic : public discr_basic
{
protected:
int find_triplet(Element * r, Element * f, Face * harmonic_face, Storage::real v[3],
Storage::real dir, Storage::real area,
Element * t[3], Storage::real coefs[3], int ptypes[3],
int allowed = AVG_HARMONIC | AVG_NONLINEAR | AVG_REVDIST, bool f_necessery = false);
public:
discr_triplet_basic(Automatizator * aut, Mesh * m, std::string tensor_name) : discr_basic(aut,m,tensor_name) {}
virtual ~discr_triplet_basic();
};
class tpfa : public discr_basic
{
INMOST_DATA_ENUM_TYPE stncl;
Tag trans;
Tag elems;
void Build(Face * f);
public:
tpfa(Automatizator * aut, Mesh *m, std::string tensor);
DiscrType GetDiscrType() const { return TPFA; }
SchemeType GetSchemeType() const { return Linear; }
expr Grad(const expr & param) const;
expr BoundaryGrad(const expr & param) const {return Grad(param);}
expr Interp(ElementType etype,const expr & param) const;
void Export(std::ostream & out, Storage::real trans_scale, Storage::real vol_scale) const;
void Init();
void Update();
bool needBuild(Face * f) const { return discr_basic::needBuild(f); }
bool needUpdate(Face * f) const { return discr_basic::needUpdate(f); }
~tpfa();
};
class tpfa2 : public discr_basic
{
INMOST_DATA_ENUM_TYPE stncl;
Tag trans;
Tag elems;
void Build(Face * f);
public:
tpfa2(Automatizator * aut, Mesh *m, std::string tensor);
DiscrType GetDiscrType() const { return TPFA; }
SchemeType GetSchemeType() const { return Linear; }
expr Grad(const expr & param) const;
expr BoundaryGrad(const expr & param) const {return Grad(param);}
expr Interp(ElementType etype,const expr & param) const;
void Export(std::ostream & out, Storage::real trans_scale, Storage::real vol_scale) const;
void Init();
void Update();
bool needBuild(Face * f) const { return discr_basic::needBuild(f); }
bool needUpdate(Face * f) const { return discr_basic::needUpdate(f); }
~tpfa2();
};
class mpfa_o : public discr_basic
{
INMOST_DATA_ENUM_TYPE stncl;
Tag trans;
Tag elems;
void Build(Face * f);
public:
mpfa_o(Automatizator * out, Mesh *m, std::string tensor);
DiscrType GetDiscrType() const { return MPFA; }
SchemeType GetSchemeType() const { return Linear; }
expr Grad(const expr & param) const;
expr BoundaryGrad(const expr & param) const {return Grad(param);}
expr Interp(ElementType etype,const expr & param) const;
void Export(std::ostream & out, Storage::real trans_scale, Storage::real vol_scale) const;
void Init();
void Update();
bool needBuild(Face * f) const { return discr_basic::needBuild(f); } //TODO
bool needUpdate(Face * f) const { return discr_basic::needBuild(f); } //TODO
~mpfa_o();
};
class basic_a : public discr_triplet_basic
{
MIDType mark_next;
INMOST_DATA_ENUM_TYPE dwn, upw, intrp, add;
INMOST_DATA_ENUM_TYPE trpr;
INMOST_DATA_ENUM_TYPE trpl;
Tag trpl_elems, trpr_elems;
Tag trpl_coefs, trpr_coefs;
Tag trans;
Tag intrp_stncl, intrp_coefs;
Tag pos;
void MarkRecursive(Element * f);
void Build(Face * f);
static void get_r(Storage * e, Automatizator::stencil_pairs & out, void * user_data);
static void get_l(Storage * e, Automatizator::stencil_pairs & out, void * user_data);
public:
basic_a(Automatizator * out, Mesh *m, std::string tensor);
void Init();
void Update();
bool needBuild(Face * f) const { return discr_basic::needBuild(f); } //TODO
bool needUpdate(Face * f) const { return discr_basic::needBuild(f); } //TODO
virtual ~basic_a();
};
class ntpfa_a : public discr_triplet_basic
{
MIDType mark_next;
INMOST_DATA_ENUM_TYPE dwn, upw, intrp;
INMOST_DATA_ENUM_TYPE trpr, addr;
INMOST_DATA_ENUM_TYPE trpl, addl;
Tag trpl_elems, trpr_elems;
Tag trpl_coefs, trpr_coefs;
Tag trpl_add, trpr_add;
Tag intrp_stncl, intrp_coefs;
Tag pos;
void MarkRecursive(Element * f);
void MarkRecursiveEdge(Element * f);
void Build(Face * f);
static void get_r(Storage * e, Automatizator::stencil_pairs & out, void * user_data);
static void get_l(Storage * e, Automatizator::stencil_pairs & out, void * user_data);
public:
ntpfa_a(Automatizator * out, Mesh *m, std::string tensor);
DiscrType GetDiscrType() const { return NTPFA; }
SchemeType GetSchemeType() const { return Nonlinear; }
expr Grad(const expr & param) const;
expr BoundaryGrad(const expr & param) const {return Grad(param);}
expr Interp(ElementType etype,const expr & param) const;
void Export(std::ostream & out, Storage::real trans_scale, Storage::real vol_scale) const;
void Init();
void Update();
bool needBuild(Face * f) const { return discr_basic::needBuild(f); } //TODO
bool needUpdate(Face * f) const { return discr_basic::needBuild(f); } //TODO
~ntpfa_a();
};
class nmpfa_a : public basic_a
{
public:
nmpfa_a(Automatizator * out, Mesh *m, std::string tensor);
DiscrType GetDiscrType() const { return NMPFA; }
SchemeType GetSchemeType() const { return Nonlinear; }
expr Grad(const expr & param) const;
expr BoundaryGrad(const expr & param) const {return Grad(param);}
expr Interp(ElementType etype,const expr & param) const;
void Export(std::ostream & out, Storage::real trans_scale, Storage::real vol_scale) const;
~nmpfa_a();
};
class ntpfa_c : public discr_triplet_basic
{
MIDType mark_next;
INMOST_DATA_ENUM_TYPE dwn, upw, intrp;
INMOST_DATA_ENUM_TYPE trpr1, trpl1, addl1, addr1;
INMOST_DATA_ENUM_TYPE trpr2, trpl2, addl2, addr2;
Tag trpl1_elems, trpl2_elems, trpr1_elems, trpr2_elems;
Tag trpl1_coefs, trpl2_coefs, trpr1_coefs, trpr2_coefs;
Tag trpl1_add, trpl2_add, trpr1_add, trpr2_add;
Tag intrp_stncl, intrp_coefs;
void MarkRecursive(Element * f);
void MarkRecursiveEdge(Element * f);
void Build(Face * f);
static void get_r(Storage * e, Automatizator::stencil_pairs & out, void * user_data);
static void get_l(Storage * e, Automatizator::stencil_pairs & out, void * user_data);
public:
ntpfa_c(Automatizator * out, Mesh *m, std::string tensor);
DiscrType GetDiscrType() const { return NTPFA; }
SchemeType GetSchemeType() const { return Nonlinear; }
expr Grad(const expr & param) const;
expr BoundaryGrad(const expr & param) const;
expr Interp(ElementType etype,const expr & param) const;
void Export(std::ostream & out, Storage::real trans_scale, Storage::real vol_scale) const;
void Init();
void Update();
bool needBuild(Face * f) const { return discr_basic::needBuild(f); } //TODO
bool needUpdate(Face * f) const { return discr_basic::needBuild(f); } //TODO
~ntpfa_c();
};
class ntpfa_b : public discr_triplet_basic
{
MIDType mark_next;
INMOST_DATA_ENUM_TYPE dwn, upw, intrp;
INMOST_DATA_ENUM_TYPE trpr, addr;
INMOST_DATA_ENUM_TYPE trpl, addl;
Tag trpl_elems, trpr_elems;
Tag trpl_coefs, trpr_coefs;
Tag trpl_add, trpr_add;
Tag intrp_stncl, intrp_coefs;
Tag pos;
void MarkRecursive(Element * f);
void MarkRecursiveEdge(Element * f);
void Build(Face * f);
static void get_r(Storage * e, Automatizator::stencil_pairs & out, void * user_data);
static void get_l(Storage * e, Automatizator::stencil_pairs & out, void * user_data);
public:
ntpfa_b(Automatizator * out, Mesh *m, std::string tensor);
DiscrType GetDiscrType() const { return NTPFA; }
SchemeType GetSchemeType() const { return Nonlinear; }
expr Grad(const expr & param) const;
expr BoundaryGrad(const expr & param) const {return Grad(param);}
expr Interp(ElementType etype,const expr & param) const;
void Export(std::ostream & out, Storage::real trans_scale, Storage::real vol_scale) const;
void Init();
void Update();
bool needBuild(Face * f) const { return discr_basic::needBuild(f); } //TODO
bool needUpdate(Face * f) const { return discr_basic::needBuild(f); } //TODO
~ntpfa_b();
};
class nmpfa : public discr_triplet_basic
{
MIDType mark_next;
INMOST_DATA_ENUM_TYPE dwn, upw, intrp, add;
INMOST_DATA_ENUM_TYPE trpr;
INMOST_DATA_ENUM_TYPE trpl;
Tag trpl_elems, trpr_elems;
Tag trpl_coefs, trpr_coefs;
Tag addval;
Tag intrp_stncl, intrp_coefs;
Tag pos;
void MarkRecursive(Element * f);
void Build(Face * f);
static void get_r(Storage * e, Automatizator::stencil_pairs & out, void * user_data);
static void get_l(Storage * e, Automatizator::stencil_pairs & out, void * user_data);
public:
nmpfa(Automatizator * out, Mesh *m, std::string tensor);
DiscrType GetDiscrType() const { return NMPFA; }
SchemeType GetSchemeType() const { return Nonlinear; }
expr Grad(const expr & param) const;
expr BoundaryGrad(const expr & param) const {return Interp(FACE,param);}
expr Interp(ElementType etype,const expr & param) const;
void Export(std::ostream & out, Storage::real trans_scale, Storage::real vol_scale) const;
void Init();
void Update();
bool needBuild(Face * f) const { return discr_basic::needBuild(f); } //TODO
bool needUpdate(Face * f) const { return discr_basic::needBuild(f); } //TODO
~nmpfa();
};
class mpfa_t : public discr_triplet_basic
{
MIDType mark_next;
INMOST_DATA_ENUM_TYPE scheme, interp;
Tag stncl, coefs, interp_stncl, interp_coefs;
void MarkRecursive(Element * f);
void Build(Face * f);
public:
mpfa_t(Automatizator * out, Mesh *m, std::string tensor);
DiscrType GetDiscrType() const { return NMPFA; }
SchemeType GetSchemeType() const { return Linear; }
expr Grad(const expr & param) const;
expr BoundaryGrad(const expr & param) const {return Grad(param);}
expr Interp(ElementType etype,const expr & param) const;
void Export(std::ostream & out, Storage::real trans_scale, Storage::real vol_scale) const;
void Init();
void Update();
bool needBuild(Face * f) const { return discr_basic::needBuild(f); } //TODO
bool needUpdate(Face * f) const { return discr_basic::needBuild(f); } //TODO
~mpfa_t();
};
class mpfa_y : public discr_triplet_basic
{
MIDType mark_next;
INMOST_DATA_ENUM_TYPE scheme, interp;
Tag stncl, coefs, interp_stncl, interp_coefs;
void MarkRecursive(Element * f);
void Build(Face * f);
public:
mpfa_y(Automatizator * out, Mesh *m, std::string tensor);
DiscrType GetDiscrType() const { return NMPFA; }
SchemeType GetSchemeType() const { return Linear; }
expr Grad(const expr & param) const;
expr BoundaryGrad(const expr & param) const {return Grad(param);}
expr Interp(ElementType etype,const expr & param) const;
void Export(std::ostream & out, Storage::real trans_scale, Storage::real vol_scale) const;
void Init();
void Update();
bool needBuild(Face * f) const { return discr_basic::needBuild(f); } //TODO
bool needUpdate(Face * f) const { return discr_basic::needBuild(f); } //TODO
~mpfa_y();
};
#endif
\ No newline at end of file
#include "discr.h"
#include <iomanip>
#define EPS 1e-12
typedef Storage::real vec[3];
static Storage::real lengthS(vec x) { return dot_prod(x, x); }
static Storage::real length(vec x) { return sqrt(lengthS(x)); }
static Storage::real normalize(Storage::real x[3]) { Storage::real r = length(x); if (r > 0.0) x[0] /= r, x[1] /= r, x[2] /= r; return r; }
basic_a::basic_a(Automatizator * aut, Mesh * m, std::string tensor)
:discr_triplet_basic(aut, m, tensor)
{
trpl_elems = m->CreateTag("A_TRP_L", DATA_REFERENCE, FACE, NONE, 4);
trpr_elems = m->CreateTag("A_TRP_R", DATA_REFERENCE, FACE, NONE, 4);
trpl_coefs = m->CreateTag("A_COEF_L", DATA_REAL, FACE, NONE, 4);
trpr_coefs = m->CreateTag("A_COEF_R", DATA_REAL, FACE, NONE, 4);
trans = m->CreateTag("A_TRANS", DATA_REAL, FACE, NONE, 1);
trpr = aut->RegisterStencil("a_right_triplet", trpr_elems, trpr_coefs);
trpl = aut->RegisterStencil("a_left_triplet", trpl_elems, trpl_coefs);
add = aut->RegisterStaticTag(trans);
intrp_stncl = m->CreateTag("INTRP_STNCL", DATA_REFERENCE, FACE | EDGE | NODE, NONE);
intrp_coefs = m->CreateTag("INTRP_COEFS", DATA_REAL, FACE | EDGE | NODE, NONE);
dwn = aut->RegisterStencil("back_cell", get_l);
upw = aut->RegisterStencil("front_cell", get_r);
intrp = aut->RegisterStencil("interp",intrp_stncl,intrp_coefs,add_markers);
}
void basic_a::MarkRecursive(Element * f)
{
if( f->GetElementType() == FACE && f->GetMarker(bnd_markers) )
{
Storage::reference_array lrtrp[2];
lrtrp[0] = f->ReferenceArrayDF(trpl_elems);
lrtrp[1] = f->ReferenceArrayDF(trpr_elems);
for(int i = 0; i < 2; i++)
for(Storage::reference_array::iterator it = lrtrp[i].begin(); it != lrtrp[i].end(); ++it) if( (*it) != NULL )
{
if( (*it)->GetElementType() != CELL && !(*it)->GetMarker(add_markers))
{
(*it)->SetMarker(add_markers);
(*it)->SetMarker(mark_next);
}
}
for(int i = 0; i < 2; i++)
for(Storage::reference_array::iterator it = lrtrp[i].begin(); it != lrtrp[i].end(); ++it) if( (*it) != NULL )
{
if( (*it)->GetMarker(mark_next) )
{
(*it)->RemMarker(mark_next);
MarkRecursive(*it);
}
}
}
else
{
Storage::reference_array arr = f->ReferenceArrayDV(intrp_stncl);
for(Storage::reference_array::iterator it = arr.begin(); it != arr.end(); ++it)
{
if( (*it)->GetElementType() != CELL && !(*it)->GetMarker(add_markers))
{
(*it)->SetMarker(add_markers);
(*it)->SetMarker(mark_next);
}
}
for(Storage::reference_array::iterator it = arr.begin(); it != arr.end(); ++it)
{
if( (*it)->GetMarker(mark_next) )
{
(*it)->RemMarker(mark_next);
MarkRecursive(*it);
}
}
}
}
void basic_a::Init()
{
if (!have_bnd) Boundary(0, FACE | EDGE);
for (INMOST_DATA_INTEGER_TYPE id = 0; id < m->MaxLocalID(FACE); ++id)
{
Face * f = m->FaceByLocalID(id);
if (f != NULL)
{
Build(f);
}
}
if( bnd_conds.isValid() )
for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it)
if( it->HaveData(bnd_conds) && it->GetMarker(bnd_markers) )
{
it->SetMarker(add_markers);
}
for(Mesh::iteratorElement it = m->BeginElement(EDGE | NODE); it != m->EndElement(); ++it)
{
Storage::real v[3], cnt[3], cnt0[3], l;
it->Centroid(cnt0);
Storage::real_array coefs = it->RealArray(intrp_coefs);
Storage::reference_array stncl = it->ReferenceArray(intrp_stncl);
adjacent<Element> cells = it->getAdjElements(it->GetMarker(bnd_markers) ? FACE : CELL);
Storage::real sum = 0;
for(int k = 0; k < cells.size(); ++k)
if( !it->GetMarker(bnd_markers) || cells[k].GetMarker(bnd_markers) )
{
cells[k].Centroid(cnt);
vec_diff(cnt0,cnt,v);
l = 1.0/(dot_prod(v,v) + 1.0e-9);
stncl.push_back(&cells[k]);
coefs.push_back(l);
sum += l;
}
for(int k = 0; k < coefs.size(); ++k) coefs[k] /= sum;
}
mark_next = m->CreateMarker();
for (Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it)
{
if( it->GetMarker(add_markers) )
MarkRecursive(&*it);
}
m->ReleaseMarker(mark_next);
for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it) it->SetMarker(add_markers);
m->SynchronizeMarker(UnknownMarker(),FACE|EDGE|NODE,SYNC_BIT_OR);
}
void basic_a::get_l(Storage * e, Automatizator::stencil_pairs & out, void * user_data)
{
Face * f = static_cast<Face *>(e);
Element * r = f->BackCell();
out.push_back(Automatizator::stencil_pair(r == NULL ? f : r, 1.0));
}
void basic_a::get_r(Storage * e, Automatizator::stencil_pairs & out, void * user_data)
{
Face * f = static_cast<Face *>(e);
Element * r = f->FrontCell();
out.push_back(Automatizator::stencil_pair(r == NULL ? f : r, 1.0));
}
void basic_a::Build(Face * f)
{
if (needBuild(f))
{
int ptypes[2][3];
Storage::real & T = f->RealDF(trans);
Storage::reference_array ltrp = f->ReferenceArrayDF(trpl_elems);
Storage::reference_array rtrp = f->ReferenceArrayDF(trpr_elems);
Storage::real_array lcoef = f->RealArrayDF(trpl_coefs);
Storage::real_array rcoef = f->RealArrayDF(trpr_coefs);
if( !f->GetMarker(bnd_markers) )
{
int ret1,ret2;
Storage::real xK[3],xL[3], gamma[3];