Skip to content
Snippets Groups Projects
Commit f921b237 authored by Jean-François Remacle's avatar Jean-François Remacle
Browse files

Added a parameter in GFaceCompound (the lin solver) 

started to create parallel topology in the mesh
parent f9f5a96a
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 148 additions and 44 deletions
Geo/GFaceCompound.cpp
+ 24
− 26
View file @ f921b237
......@@ -457,10 +457,24 @@ void GFaceCompound::computeALoop(std::set<GEdge*> &_unique, std::list<GEdge*> &l
GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &U0, std::list<GEdge*> &V0,
std::list<GEdge*> &U1, std::list<GEdge*> &V1)
: GFace(m, tag), _compound(compound), _U0(U0), _U1(U1), _V0(V0), _V1(V1), oct(0)
std::list<GEdge*> &U1, std::list<GEdge*> &V1,
gmshLinearSystem<double> * lsys)
: GFace(m, tag), _compound(compound), _U0(U0), _U1(U1), _V0(V0), _V1(V1), oct(0),
_lsys(lsys)
{
#ifdef HAVE_GMM
if (!_lsys) {
gmshLinearSystemGmm<double> *_lsysb = new gmshLinearSystemGmm<double>;
//gmshLinearSystemCSRGmm<double> lsys;
_lsysb->setPrec(1.e-15);
if(Msg::GetVerbosity() == 99) _lsysb->setNoisy(2);
_lsys = _lsysb;
}
#else
_lsys = new gmshLinearSystemFull<double>;
#endif
for(std::list<GFace*>::iterator it = _compound.begin(); it != _compound.end(); ++it){
if(!(*it)){
Msg::Error("Incorrect face in compound surface %d\n", tag);
......@@ -651,15 +665,7 @@ void GFaceCompound::parametrize(iterationStep step) const
{
Msg::Info("Parametrizing Surface %d coordinate (ITER %d)", tag(), step);
#ifdef HAVE_GMM
gmshLinearSystemGmm<double> lsys;
//gmshLinearSystemCSRGmm<double> lsys;
lsys.setPrec(1.e-15);
if(Msg::GetVerbosity() == 99) lsys.setNoisy(2);
#else
gmshLinearSystemFull<double> lsys;
#endif
gmshAssembler<double> myAssembler(&lsys);
gmshAssembler<double> myAssembler(_lsys);
gmshFunction<double> ONE(1.0);
......@@ -728,7 +734,7 @@ void GFaceCompound::parametrize(iterationStep step) const
}
Msg::Debug("Assembly done");
lsys.systemSolve();
_lsys->systemSolve();
Msg::Debug("System solved");
......@@ -755,8 +761,8 @@ void GFaceCompound::parametrize(iterationStep step) const
else{
itf->second[step]= value;
}
}
_lsys->clear();
}
void GFaceCompound::parametrize_conformal() const
......@@ -765,16 +771,8 @@ void GFaceCompound::parametrize_conformal() const
Msg::Info("Parametrizing Surface %d", tag());
#ifdef HAVE_GMM
gmshLinearSystemGmm<double> lsys;
//gmshLinearSystemCSRGmm<double> lsys;
lsys.setPrec(1.e-15);
lsys.setGmres(1);
if(Msg::GetVerbosity() == 99) lsys.setNoisy(2);
#else
gmshLinearSystemFull<double> lsys;
#endif
gmshAssembler<double> myAssembler(&lsys);
gmshAssembler<double> myAssembler(_lsys);
std::vector<MVertex*> ordered;
std::vector<double> coords;
bool success = orderVertices(_U0, ordered, coords);
......@@ -863,7 +861,7 @@ void GFaceCompound::parametrize_conformal() const
}
Msg::Debug("Assembly done");
lsys.systemSolve();
_lsys->systemSolve();
Msg::Debug("System solved");
it = _compound.begin();
......@@ -885,9 +883,9 @@ void GFaceCompound::parametrize_conformal() const
}
// checkOrientation();
computeNormals();
buildOct();
_lsys->clear();
}
void GFaceCompound::computeNormals(std::map<MVertex*,SVector3> &normals) const
......
Geo/GFaceCompound.h
+ 4
− 1
View file @ f921b237
......@@ -12,6 +12,7 @@
#include "GEdge.h"
#include "GEdgeCompound.h"
#include "meshGFaceOptimize.h"
#include "gmshLinearSystem.h"
/*
A GFaceCompound is a model face that is the compound of model faces.
......@@ -73,11 +74,13 @@ class GFaceCompound : public GFace {
virtual double curvature(MTriangle *t, double u, double v) const;
void printStuff() const;
bool trivial() const ;
gmshLinearSystem <double> *_lsys;
public:
typedef enum {UNITCIRCLE, SQUARE} typeOfIsomorphism;
GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &U0, std::list<GEdge*> &U1,
std::list<GEdge*> &V0, std::list<GEdge*> &V1);
std::list<GEdge*> &V0, std::list<GEdge*> &V1,
gmshLinearSystem<double>* lsys =0);
virtual ~GFaceCompound();
Range<double> parBounds(int i) const
{ return trivial() ? (*(_compound.begin()))->parBounds(i) : Range<double>(-1, 1); }
......
Geo/GeoInterpolation.cpp
+ 1
− 1
View file @ f921b237
......@@ -188,7 +188,7 @@ static void basisFuns(double u, int i, int deg, float *U, double *N)
static Vertex InterpolateNurbs(Curve *Curve, double u, int derivee)
{
static double Nb[1000];
double Nb[1000];
int span = findSpan(u, Curve->degre, List_Nbr(Curve->Control_Points), Curve->k);
basisFuns(u, span, Curve->degre, Curve->k, Nb);
Vertex p;
......
Mesh/meshPartition.cpp
+ 85
− 0
View file @ f921b237
......@@ -13,6 +13,10 @@
#include "meshPartition.h"
#include "meshPartitionObjects.h"
#include "meshPartitionOptions.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "MElementCut.h"
#include "partitionEdge.h"
//--Prototype for Chaco interface
......@@ -643,6 +647,87 @@ struct MatchBoElemToGrVertex
}
};
template <class ITERATOR>
void fillit_ ( std::multimap<MEdge,MElement*,Less_Edge> &edgeToElement,
ITERATOR it_beg, ITERATOR it_end){
for (ITERATOR IT = it_beg; IT != it_end ; ++IT){
MElement *el = *IT;
for(int j=0;j < el->getNumEdges();j++) {
MEdge e = el->getEdge(j);
edgeToElement.insert(std::make_pair(e,el));
}// for every edge of the elem
}// for every elem of the face
}
template <class ITERATOR>
void fillit_ ( std::multimap<MVertex*,MElement*> &vertexToElement,
ITERATOR it_beg, ITERATOR it_end){
for (ITERATOR IT = it_beg; IT != it_end ; ++IT){
MElement *el = *IT;
for(int j=0;j < el->getNumVertices();j++) {
MVertex* e = el->getVertex(j);
vertexToElement.insert(std::make_pair(e,el));
}// for every edge of the elem
}// for every elem of the face
}
void assignPartitionBoundary (GModel *model,
MEdge &me,
std::set<partitionEdge*, Less_partitionEdge> &pedges,
std::vector<MElement*> &v){
std::vector<int> v2;
v2.push_back(v[0]->getPartition());
for (int i=1;i<v.size();i++){
bool found = false;
for (int j=0;j<v2.size();j++){
if (v[i]->getPartition() == v2[j]){
found = true;
break;
}
}
if (!found)v2.push_back(v[i]->getPartition());
}
if (v2.size() < 2)return;
partitionEdge pe (model, 1, 0, 0, v2);
std::set<partitionEdge*, Less_partitionEdge>::iterator it = pedges.find(&pe);
partitionEdge *ppe;
if (it == pedges.end()){
ppe = new partitionEdge(model, -pedges.size()-1, 0, 0, v2);
pedges.insert (ppe);
model->add(ppe);
}
else ppe = *it;
ppe->lines.push_back(new MLine (me.getVertex(0),me.getVertex(1)));
}
int CreatePartitionBoundaries (GModel *model) {
std::multimap<MEdge,MElement*,Less_Edge> edgeToElement;
std::set<partitionEdge*, Less_partitionEdge> pedges;
for(GModel::fiter it = model->firstFace(); it != model->lastFace(); ++it){
fillit_ ( edgeToElement,(*it)->triangles.begin(),(*it)->triangles.end());
fillit_ ( edgeToElement,(*it)->quadrangles.begin(),(*it)->quadrangles.end());
fillit_ ( edgeToElement,(*it)->polygons.begin(),(*it)->polygons.end());
}
{
std::multimap<MEdge,MElement*,Less_Edge>::iterator it = edgeToElement.begin();
Equal_Edge oper;
while ( it != edgeToElement.end() ){
MEdge e = it->first;
std::vector<MElement*> voe;
do {
voe.push_back(it->second);
++it;
}while (! oper (e,it->first));
assignPartitionBoundary (model,e,pedges,voe);
}
}
}
/*******************************************************************************
*
......
Mesh/meshPartition.h
+ 1
− 0
View file @ f921b237
......@@ -25,5 +25,6 @@ int MakeGraph(GModel *const model, Graph &graph,
BoElemGrVec *const boElemGrVec = 0);
int PartitionGraph(Graph &graph, meshPartitionOptions &options);
int PartitionMesh(GModel *const model, meshPartitionOptions &options);
int CreatePartitionBoundaries (GModel *model);
#endif
Numeric/gmshLinearSystem.h
+ 1
− 0
View file @ f921b237
......@@ -16,6 +16,7 @@ class gmshLinearSystem {
virtual ~gmshLinearSystem (){}
virtual bool isAllocated() const = 0;
virtual void allocate(int nbRows) = 0;
virtual void clear() = 0;
virtual void addToMatrix(int _row, int _col, scalar val) = 0;
virtual scalar getFromMatrix(int _row, int _col) const = 0;
virtual void addToRightHandSide(int _row, scalar val) = 0;
......
Numeric/gmshLinearSystemCSR.h
+ 4
− 0
View file @ f921b237
......@@ -37,6 +37,10 @@ class gmshLinearSystemCSR : public gmshLinearSystem<scalar> {
: sorted(false), a_(0) {}
virtual bool isAllocated() const { return a_ != 0; }
virtual void allocate(int) ;
virtual void clear()
{
allocate(0);
}
virtual ~gmshLinearSystemCSR()
{
allocate(0);
......
Numeric/gmshLinearSystemFull.h
+ 13
− 6
View file @ f921b237
......@@ -24,19 +24,26 @@ class gmshLinearSystemFull : public gmshLinearSystem<scalar> {
virtual bool isAllocated() const { return _a != 0; }
virtual void allocate(int _nbRows)
{
if(_a) delete _a;
if(_x) delete _x;
if(_b) delete _b;
clear();
_a = new gmshMatrix<scalar>(_nbRows, _nbRows);
_b = new gmshVector<scalar>(_nbRows);
_x = new gmshVector<scalar>(_nbRows);
}
virtual ~gmshLinearSystemFull()
{
clear();
}
virtual void clear()
{
if(_a){
delete _a;
delete _b;
delete _x;
}
_a = 0;
}
virtual void addToMatrix(int _row, int _col, scalar _val)
{
if(_val != 0.0) (*_a)(_row, _col) += _val;
......
Numeric/gmshLinearSystemGmm.h
+ 12
− 6
View file @ f921b237
......@@ -29,19 +29,25 @@ class gmshLinearSystemGmm : public gmshLinearSystem<scalar> {
virtual bool isAllocated() const { return _a != 0; }
virtual void allocate(int _nbRows)
{
if(_a) delete _a;
if(_x) delete _x;
if(_b) delete _b;
clear();
_a = new gmm::row_matrix< gmm::wsvector<scalar> >(_nbRows, _nbRows);
_b = new std::vector<scalar>(_nbRows);
_x = new std::vector<scalar>(_nbRows);
}
virtual ~gmshLinearSystemGmm()
{
clear();
}
virtual void clear()
{
if (_a){
delete _a;
delete _b;
delete _x;
}
_a = 0;
}
virtual void addToMatrix(int _row, int _col, scalar _val)
{
if(_val != 0.0) (*_a)(_row, _col) += _val;
......
benchmarks/3d/Torus.geo
+ 1
− 1
View file @ f921b237
......@@ -13,4 +13,4 @@ Line Loop(5) = {4,1,2,3};
Plane Surface(6) = {5};
Extrude Surface{6, {0.0,1,0}, {0,0.0,0.0}, 1*3.14159/2};
Recombine Surface {6, 27, 23, 15, 19, 28};
//Recombine Surface {6, 27, 23, 15, 19, 28};
benchmarks/boolean/FILLET.geo
+ 2
− 3
View file @ f921b237
......@@ -3,7 +3,9 @@ H = 30;
Z = 10;
OCCShape("Box",{0,0,0,L,H,Z},"none");
R = 4;
OCCShape("Fillet",{1:12,R},"none");
/*
OCCShape("Cone",{0*L/2,H/2,-Z,0,0,1,.3*R,2*R,3*Z},"Fuse");
OCCShape("Fillet",{1,R},"none");
OCCShape("Fillet",{14,R/8},"none");
......@@ -11,6 +13,3 @@ OCCShape("Cone",{0.99*L/2,H/2,-Z,0,0,1,.3*R,2*R,3*Z},"Fuse");
OCCShape("Fillet",{1,R},"none");
OCCShape("Fillet",{83,R/8},"none");
OCCShape("Cone",{0.99*L,H/2,-Z,0,0,1,3*R,.2*R,3*Z},"Cut");
// not ready yet for seams!
//Compound Surface(100000000) = {31, 40, 43, 13};
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment