From 511e46255edceba7f233b2dc018cf335b9e74768 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Wed, 6 Sep 2006 18:42:24 +0000
Subject: [PATCH] - start work on fourierModel - add support for pyramids in
BDF (it's not in the official spec I have here, but all the files I found
on the web use CPYRAM)
---
Geo/GModelIO.cpp | 10 ++-
Geo/MElement.h | 2 +-
Geo/MVertex.h | 19 +++++-
Geo/fourierModel.cpp | 155 +++++++++++++++++++++++++++++++++++++------
4 files changed, 162 insertions(+), 24 deletions(-)
diff --git a/Geo/GModelIO.cpp b/Geo/GModelIO.cpp
index e0b64aa7b2..bb1d95e189 100644
--- a/Geo/GModelIO.cpp
+++ b/Geo/GModelIO.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO.cpp,v 1.42 2006-09-04 15:58:22 geuzaine Exp $
+// $Id: GModelIO.cpp,v 1.43 2006-09-06 18:42:24 geuzaine Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -1590,7 +1590,7 @@ int GModel::readBDF(const std::string &name)
char buffer[256];
std::map<int, MVertex*> vertices;
- std::map<int, std::vector<MElement*> > elements[6];
+ std::map<int, std::vector<MElement*> > elements[7];
// nodes can be defined after elements, so parse the file twice
@@ -1643,6 +1643,10 @@ int GModel::readBDF(const std::string &name)
if(readElementBDF(fp, buffer, 6, 6, &num, ®ion, n, vertices))
elements[5][region].push_back(new MPrism(n, num));
}
+ else if(!strncmp(buffer, "CPYRAM", 6)){
+ if(readElementBDF(fp, buffer, 6, 5, &num, ®ion, n, vertices))
+ elements[6][region].push_back(new MPyramid(n, num));
+ }
}
}
@@ -1698,6 +1702,8 @@ int GModel::writeBDF(const std::string &name, int format, bool saveAll,
(*it)->hexahedra[i]->writeBDF(fp, format, (*it)->tag());
for(unsigned int i = 0; i < (*it)->prisms.size(); i++)
(*it)->prisms[i]->writeBDF(fp, format, (*it)->tag());
+ for(unsigned int i = 0; i < (*it)->pyramids.size(); i++)
+ (*it)->pyramids[i]->writeBDF(fp, format, (*it)->tag());
}
fprintf(fp, "ENDDATA\n");
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 7f639644ab..daba42afc1 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -717,7 +717,7 @@ class MPyramid : public MElement {
virtual int getTypeForMSH(){ return PYR1; }
virtual int getTypeForUNV(){ return 0; } // not available
virtual char *getStringForPOS(){ return "SY"; }
- virtual char *getStringForBDF(){ return 0; } // not available
+ virtual char *getStringForBDF(){ return "CPYRAM"; }
virtual int getVolumeSign()
{
double mat[3][3];
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index c6f32f211a..101430ed72 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -76,6 +76,9 @@ class MVertex{
// Get ith parameter
virtual bool getParameter(int i, double &par){ return false; }
+ // Get the data associated with this vertex
+ virtual void *getData(){ return 0; }
+
// IO routines
void writeMSH(FILE *fp, bool binary=false, double scalingFactor=1.0);
void writeMSH(FILE *fp, double version, bool binary, int num,
@@ -99,7 +102,7 @@ class MEdgeVertex : public MVertex{
};
class MFaceVertex : public MVertex{
- private:
+ protected:
double _u, _v;
public :
MFaceVertex(double x, double y, double z, GEntity *ge, double u, double v)
@@ -110,6 +113,20 @@ class MFaceVertex : public MVertex{
virtual bool getParameter(int i, double &par){ par = (i ? _v : _u); return true; }
};
+template<class T>
+class MDataFaceVertex : public MFaceVertex{
+ private:
+ T _data;
+ public :
+ MDataFaceVertex(double x, double y, double z, GEntity *ge, double u, double v, T data)
+ : MFaceVertex(x, y, z, ge, u, v), _data(data)
+ {
+ }
+ ~MDataFaceVertex(){}
+ virtual bool getData(T &data){ data = _data; return true; }
+ virtual void *getData(){ return (void*)&_data; }
+};
+
class MVertexLessThanLexicographic{
public:
static double tolerance;
diff --git a/Geo/fourierModel.cpp b/Geo/fourierModel.cpp
index 3bba29d23e..1667367664 100644
--- a/Geo/fourierModel.cpp
+++ b/Geo/fourierModel.cpp
@@ -1,5 +1,9 @@
#include "fourierModel.h"
#include "Message.h"
+#include "Context.h"
+#include "Views.h"
+
+extern Context_T CTX;
#if defined(HAVE_FOURIER_MODEL)
@@ -7,6 +11,111 @@
static model *FM = 0;
+class meshFourierFace{
+public:
+ void operator() (GFace *gf)
+ {
+ int M = 30, N = 30;
+ for(int i = 0; i < M; i++){
+ for(int j = 0; j < N; j++){
+ double u = i/(double)(M - 1);
+ double v = j/(double)(N - 1);
+ GPoint p = gf->point(u, v);
+ double pou;
+ FM->GetPou(gf->tag(), u, v, pou);
+ gf->mesh_vertices.push_back
+ (new MDataFaceVertex<double>(p.x(), p.y(), p.z(), gf, u, v, pou));
+ }
+ }
+ for(int i = 0; i < M - 1; i++){
+ for(int j = 0; j < N - 1; j++){
+ gf->quadrangles.push_back(new MQuadrangle(gf->mesh_vertices[i * N + j],
+ gf->mesh_vertices[i * N + (j + 1)],
+ gf->mesh_vertices[(i + 1) * N + (j + 1)],
+ gf->mesh_vertices[(i + 1) * N + j]));
+ }
+ }
+ }
+};
+
+
+static void getParamForPointInOverlaps(GModel *m, std::vector<int> &overlaps,
+ double x, double y, double z,
+ std::vector<std::pair<GFace*, SPoint2> > ¶m)
+{
+ // we only loop on patches that are known a priori to overlap with
+ // the current patch (this way we don't normalize the POU across
+ // hard edges)
+ const double tol = 1.e-2; // need to adapt this w.r.t size of model
+ for(unsigned int i = 0; i < overlaps.size(); i++){
+ GFace *f = m->faceByTag(overlaps[i]);
+ SPoint2 uv = f->parFromPoint(SPoint3(x, y, z));
+ if(f->containsParam(uv)){
+ GPoint xyz = f->point(uv);
+ if(fabs(xyz.x() - x) < tol && fabs(xyz.y() - y) < tol && fabs(xyz.z() - z) < tol)
+ param.push_back(std::pair<GFace*, SPoint2>(f, uv));
+ }
+ }
+}
+
+class computePartitionOfUnity{
+public:
+ void operator() (GFace *gf)
+ {
+ std::vector<int> overlaps;
+ FM->GetNeighbor(gf->tag(), overlaps);
+ for(unsigned int i = 0; i < gf->mesh_vertices.size(); i++){
+ MVertex *v = gf->mesh_vertices[i];
+ std::vector<std::pair<GFace*, SPoint2> > param;
+ getParamForPointInOverlaps(gf->model(), overlaps, v->x(), v->y(), v->z(), param);
+ if(param.empty()){
+ Msg(GERROR, "Point (%g,%g,%g) does not belong to any patch",
+ v->x(), v->y(), v->z());
+ }
+ else{
+ double allPou = 0.;
+ for(unsigned int i = 0; i < param.size(); i++){
+ double u2 = param[i].second[0], v2 = param[i].second[1];
+ double pou2;
+ FM->GetPou(param[i].first->tag(), u2, v2, pou2);
+ allPou += pou2;
+ }
+ double *pou = (double*)v->getData();
+ *pou = *pou / allPou;
+ }
+ }
+ }
+};
+
+class exportFourierFace{
+public:
+ void operator() (GFace *gf)
+ {
+ Post_View *view = BeginView(1);
+ for(unsigned int i = 0; i < gf->quadrangles.size(); i++){
+ MElement *e = gf->quadrangles[i];
+ double x[4], y[4], z[4], val[4];
+ for(int j = 0; j < 4; j++){
+ MVertex *v = e->getVertex(j);
+ x[j] = v->x();
+ y[j] = v->y();
+ z[j] = v->z();
+ val[j] = *(double*)v->getData();
+ }
+ for(int j = 0; j < 4; j++) List_Add(view->SQ, &x[j]);
+ for(int j = 0; j < 4; j++) List_Add(view->SQ, &y[j]);
+ for(int j = 0; j < 4; j++) List_Add(view->SQ, &z[j]);
+ for(int j = 0; j < 4; j++) List_Add(view->SQ, &val[j]);
+ view->NbSQ++;
+ }
+ char name[1024], filename[1024];
+ sprintf(name, "patch%d", gf->tag());
+ sprintf(filename, "patch%d", gf->tag());
+ EndView(view, 1, filename, name);
+ }
+};
+
+
fourierModel::fourierModel(const std::string &name)
: GModel(name)
{
@@ -14,10 +123,21 @@ fourierModel::fourierModel(const std::string &name)
Msg(INFO, "Fourier model created: %d patches", FM->GetNumPatches());
+ // create one face per patch
for(int i = 0; i < FM->GetNumPatches(); i++)
add(new fourierFace(this, i));
+
+ // mesh each face
+ std::for_each(firstFace(), lastFace(), meshFourierFace());
+
+ // compute partition of unity
+ std::for_each(firstFace(), lastFace(), computePartitionOfUnity());
+
+ // visualize as a post-pro view
+ std::for_each(firstFace(), lastFace(), exportFourierFace());
}
+
fourierModel::~fourierModel()
{
delete FM;
@@ -27,26 +147,13 @@ fourierModel::~fourierModel()
fourierFace::fourierFace(GModel *m, int num)
: GFace(m, num)
{
- int M = 30, N = 30;
- for(int i = 0; i < M; i++){
- for(int j = 0; j < N; j++){
- GPoint p = point(i/(double)(M - 1), j/(double)(N - 1));
- mesh_vertices.push_back(new MVertex(p.x(), p.y(), p.z(), this));
- }
- }
- for(int i = 0; i < M - 1; i++){
- for(int j = 0; j < N - 1; j++){
- quadrangles.push_back(new MQuadrangle(mesh_vertices[i * M + j],
- mesh_vertices[i * M + (j + 1)],
- mesh_vertices[(i + 1) * M + (j + 1)],
- mesh_vertices[(i + 1) * M + j]));
- }
- }
}
Range<double> fourierFace::parBounds(int i) const
{
- return Range<double>(0., 1.);
+ double min, max;
+ FM->GetParamBounds(tag(), i, min, max);
+ return Range<double>(min, max);
}
GPoint fourierFace::point(double par1, double par2) const
@@ -58,7 +165,7 @@ GPoint fourierFace::point(double par1, double par2) const
GPoint fourierFace::point(const SPoint2 &pt) const
{
- throw;
+ return point(pt[0], pt[1]);
}
GPoint fourierFace::closestPoint(const SPoint3 & queryPoint)
@@ -73,7 +180,13 @@ int fourierFace::containsPoint(const SPoint3 &pt) const
int fourierFace::containsParam(const SPoint2 &pt) const
{
- throw;
+ double umin, umax, vmin, vmax;
+ FM->GetParamBounds(tag(), 0, umin, umax);
+ FM->GetParamBounds(tag(), 1, vmin, vmax);
+ const double tol = 1.e-6;
+ if(pt[0] < umin - tol || pt[0] > umax + tol) return 0;
+ if(pt[1] < vmin - tol || pt[1] > vmax + tol) return 0;
+ return 1;
}
SVector3 fourierFace::normal(const SPoint2 ¶m) const
@@ -86,9 +199,11 @@ GEntity::GeomType fourierFace::geomType() const
return GEntity::DiscreteSurface;
}
-SPoint2 fourierFace::parFromPoint(const SPoint3 &) const
+SPoint2 fourierFace::parFromPoint(const SPoint3 &p) const
{
- throw;
+ double u, v;
+ FM->GetParamFromPoint(tag(), p.x(), p.y(), p.z(), u, v);
+ return SPoint2(u, v);
}
#endif
--
GitLab