diff --git a/Common/Options.cpp b/Common/Options.cpp
index c96f6face5f37f7f80a83594db691f30f31047b0..8b597cb5654c66893b331c1aae83373c5e67f77f 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -5176,12 +5176,10 @@ double opt_mesh_recombine_all(OPT_ARGS_NUM)
double opt_mesh_recombine3d_all(OPT_ARGS_NUM)
{
-#if defined(HAVE_FLTK)
if(action & GMSH_SET){
CTX::instance()->mesh.recombine3DAll = (int)val;
-#if defined(HAVE_FLTK)
-#endif
}
+#if defined(HAVE_FLTK)
if(FlGui::available() && (action & GMSH_GUI)){
FlGui::instance()->options->mesh.butt[22]->value
(CTX::instance()->mesh.recombine3DAll);
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 40c7d66ba97d6ad95c043621e3b022a9c91719b1..e24a6ff8be6f3cf12d8778751f589469595eeba8 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -93,6 +93,11 @@ class GFace : public GEntity
// edges that bound the face
virtual std::list<GEdge*> edges() const { return l_edges; }
virtual std::list<int> edgeOrientations() const { return l_dirs; }
+ inline bool containsEdge (int iEdge) const {
+ for (std::list<GEdge*>::const_iterator it = l_edges.begin() ; it !=l_edges.end() ; ++it)
+ if ((*it)->tag() == iEdge) return true;
+ return false;
+ }
// edges that are embedded in the face
virtual std::list<GEdge*> embeddedEdges() const { return embedded_edges; }
diff --git a/Geo/GModelIO_MESH.cpp b/Geo/GModelIO_MESH.cpp
index f15e22d4b02a2c62cd0049076ad2b45232035d12..7f43d9c57c57bffd32434baec211aca3448e1f60 100644
--- a/Geo/GModelIO_MESH.cpp
+++ b/Geo/GModelIO_MESH.cpp
@@ -109,8 +109,7 @@ int GModel::readMESH(const std::string &name)
for(int i = 0; i < nbe; i++) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
int n[6], cl;
- sscanf(buffer, "%d %d %d %d %d %d %d", &n[0], &n[1], &n[2],
- &n[3], &n[4], &n[5], &cl);
+ sscanf(buffer, "%d %d %d %d %d %d %d", &n[0], &n[1], &n[2], &n[3], &n[4], &n[5],&cl);
for(int j = 0; j < 6; j++) n[j]--;
std::vector<MVertex*> vertices;
if(!getVertices(6, n, vertexVector, vertices)) return 0;
@@ -155,8 +154,7 @@ int GModel::readMESH(const std::string &name)
for(int i = 0; i < nbe; i++) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
int n[10], cl;
- sscanf(buffer, "%d %d %d %d %d %d %d %d %d %d %d", &n[0], &n[1], &n[2],
- &n[3], &n[4], &n[5], &n[6], &n[7], &n[9], &n[8], &cl);
+ sscanf(buffer, "%d %d %d %d %d %d %d %d %d %d %d", &n[0], &n[1], &n[2], &n[3],&n[4], &n[5], &n[6], &n[7], &n[9], &n[8], &cl);
for(int j = 0; j < 10; j++) n[j]--;
std::vector<MVertex*> vertices;
if(!getVertices(10, n, vertexVector, vertices)) return 0;
@@ -204,7 +202,7 @@ int GModel::writeMESH(const std::string &name, int elementTagType,
int numVertices = indexMeshVertices(saveAll);
- fprintf(fp, " MeshVersionFormatted 1\n");
+ fprintf(fp, " MeshVersionFormatted 2\n");
fprintf(fp, " Dimension\n");
fprintf(fp, " 3\n");
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index cb40e80db9f39cc868abeda76984043f657d1a29..8b30928fa84ef7af39228e2da13e9fc832ffae09 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -631,6 +631,7 @@ Volume *Create_Volume(int Num, int Typ)
pV->SurfacesOrientations = List_Create(1, 2, sizeof(int));
pV->SurfacesByTag = List_Create(1, 2, sizeof(int));
pV->Extrude = NULL;
+ pV->Recombine3D = 0;
return pV;
}
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 9a5f07bfe76cd206f99c747942a0ae075065c82a..0ccec167fbcf3b871adf0ec771db8927279afed5 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1067,7 +1067,12 @@ void MElement::writeMESH(FILE *fp, int elementTagType, int elementary,
{
setVolumePositive();
for(int i = 0; i < getNumVertices(); i++)
- fprintf(fp, " %d", getVertex(i)->getIndex());
+ if (getTypeForMSH() == MSH_TET_10 && i == 8)
+ fprintf(fp, " %d", getVertex(9)->getIndex());
+ else if (getTypeForMSH() == MSH_TET_10 && i == 9)
+ fprintf(fp, " %d", getVertex(8)->getIndex());
+ else
+ fprintf(fp, " %d", getVertex(i)->getIndex());
fprintf(fp, " %d\n", (elementTagType == 3) ? _partition :
(elementTagType == 2) ? physical : elementary);
}
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index 22b32a04d6f5dd59fac3c1abcd017de249e9b8a7..ea18b06527c3d73fb9d71b6af0f483b7d2c6421c 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -15,6 +15,7 @@
#include "Context.h"
#include "MTriangle.h"
#include "VertexArray.h"
+#include "Field.h"
gmshFace::gmshFace(GModel *m, Surface *face)
: GFace(m, face->Num), s(face)
@@ -330,7 +331,7 @@ bool gmshFace::containsPoint(const SPoint3 &pt) const
bool gmshFace::buildSTLTriangulation(bool force)
{
return false;
- /*
+
if(va_geom_triangles){
if(force)
delete va_geom_triangles;
@@ -341,31 +342,34 @@ bool gmshFace::buildSTLTriangulation(bool force)
stl_vertices.clear();
stl_triangles.clear();
- contextMeshOptions _temp = CTX::instance()->mesh;
- CTX::instance()->mesh.lcFromPoints = 0;
- CTX::instance()->mesh.lcFromCurvature = 1;
- CTX::instance()->mesh.lcExtendFromBoundary = 0;
- CTX::instance()->mesh.scalingFactor = 1;
- CTX::instance()->mesh.lcFactor = 1;
- CTX::instance()->mesh.order = 1;
- CTX::instance()->mesh.lcIntegrationPrecision = 1.e-5;
- std::for_each(l_edges.begin(), l_edges.end(), meshGEdge(true));
- meshGFace mesher (false,true);
- mesher(this);
- printf("%d triangles face %d\n",triangles_stl.size(),tag());
- CTX::instance()->mesh = _temp;
+ if (!triangles.size()){
+ contextMeshOptions _temp = CTX::instance()->mesh;
+ FieldManager *fields = model()->getFields();
+ int BGM = fields->getBackgroundField();
+ fields->setBackgroundField(0);
+ CTX::instance()->mesh.lcFromPoints = 0;
+ CTX::instance()->mesh.lcFromCurvature = 1;
+ CTX::instance()->mesh.lcExtendFromBoundary = 0;
+ CTX::instance()->mesh.scalingFactor = 1;
+ CTX::instance()->mesh.lcFactor = 1;
+ CTX::instance()->mesh.order = 1;
+ CTX::instance()->mesh.lcIntegrationPrecision = 1.e-3;
+ // CTX::instance()->mesh.Algorithm = 5;
+ model()->mesh(2);
+ CTX::instance()->mesh = _temp;
+ fields->setBackgroundField(fields->get(BGM));
+ }
std::map<MVertex*,int> _v;
int COUNT =0;
- for (int j = 0; j < triangles_stl.size(); j++){
- int C[3];
+ for (int j = 0; j < triangles.size(); j++){
for (int i=0;i<3;i++){
std::map<MVertex*,int>::iterator it =
- _v.find(triangles_stl[j]->getVertex(j));
+ _v.find(triangles[j]->getVertex(j));
if (it != _v.end()){
stl_triangles.push_back(COUNT);
- _v[triangles_stl[j]->getVertex(j)] = COUNT++;
+ _v[triangles[j]->getVertex(j)] = COUNT++;
}
else stl_triangles.push_back(it->second);
}
@@ -396,5 +400,5 @@ bool gmshFace::buildSTLTriangulation(bool force)
}
va_geom_triangles->finalize();
return true;
- */
+
}
diff --git a/Geo/gmshRegion.cpp b/Geo/gmshRegion.cpp
index fcec42f11be6a8aebe08b96dd90959cfbcdc8ed7..499fcca1f3ee40ec8e96ab448f522c6b975d5bd0 100644
--- a/Geo/gmshRegion.cpp
+++ b/Geo/gmshRegion.cpp
@@ -39,6 +39,7 @@ gmshRegion::gmshRegion(GModel *m, ::Volume *volume)
Msg::Error("Unknown surface %d", is);
}
resetMeshAttributes();
+
}
void gmshRegion::resetMeshAttributes()
diff --git a/Graphics/drawPost.cpp b/Graphics/drawPost.cpp
index 08d19d862d8c59b1f068395ee89be836d09320f9..0e2f9671e0a4a6d00b5143c109dbacc8695c9f30 100644
--- a/Graphics/drawPost.cpp
+++ b/Graphics/drawPost.cpp
@@ -68,8 +68,19 @@ static void drawArrays(drawContext *ctx, PView *p, VertexArray *va, GLint type,
double v0 = char2float(*n0), v1 = char2float(*n1);
ctx->drawTaperedCylinder(opt->lineWidth, v0, v1, 0., 1., x, y, z, opt->light);
}
- else
+ else if (opt->lineType == 1)
ctx->drawCylinder(opt->lineWidth, x, y, z, opt->light);
+ else { /// 2D (for now) MNT DIAGRAMS FOR FRAMES
+ float l = sqrt ((p0[0]-p1[0])*(p0[0]-p1[0]) + (p0[1]-p1[1])*(p0[1]-p1[1]) + (p0[2]-p1[2])*(p0[2]-p1[2]) );
+ char *n0 = va->getNormalArray(3 * i);
+ char *n1 = va->getNormalArray(3 * (i + 1));
+ double v0 = char2float(*n0), v1 = char2float(*n1);
+ float dir [3] = {(p1[0]-p0[0])/l , (p1[1]-p0[1])/l , (p1[2]-p0[2])/l };
+ printf("%g %g %g %g %g %g\n",v0,v1,p0[0],p0[1],p1[0],p1[1]);
+ ctx->drawVector(1, 0,
+ p0[0] - dir[1] * v0 , p0[1] + dir[0] * v0 , 0.0,
+ p1[0] - dir[1] * v1 , p1[1] + dir[0] * v1 , 0.0,opt->light);
+ }
}
}
else{
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index 20289b0b80fe0a9cf4bc7de12ea206cd877aa398..2c4ca6d83174f4d2a4ca11a8c61cbc744ab87986 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -1667,20 +1667,16 @@ class AttractorField : public Field
offset.push_back(0);
for(std::list<int>::iterator it = faces_id.begin();
it != faces_id.end(); ++it) {
- Surface *s = FindSurface(*it);
- if(!s) {
- GFace *f = GModel::current()->getFaceByTag(*it);
- if (f){
- SBoundingBox3d bb = f->bounds();
- SVector3 dd = bb.max() - bb.min();
- double maxDist = dd.norm() / n_nodes_by_edge ;
- f->fillPointCloud(maxDist, &points, &uvpoints);
- offset.push_back(points.size());
- }
+ GFace *f = GModel::current()->getFaceByTag(*it);
+ if (f){
+ SBoundingBox3d bb = f->bounds();
+ SVector3 dd = bb.max() - bb.min();
+ double maxDist = dd.norm() / n_nodes_by_edge ;
+ f->fillPointCloud(maxDist, &points, &uvpoints);
+ offset.push_back(points.size());
}
}
-
-
+
int totpoints =
nodes_id.size() +
(n_nodes_by_edge-2) * edges_id.size() +
@@ -1697,47 +1693,27 @@ class AttractorField : public Field
int k = 0;
for(std::list<int>::iterator it = nodes_id.begin();
it != nodes_id.end(); ++it) {
- Vertex *v = FindPoint(*it);
- if(v) {
- getCoord(v->Pos.X, v->Pos.Y, v->Pos.Z, zeronodes[k][0],
- zeronodes[k][1], zeronodes[k][2]);
+ GVertex *gv = GModel::current()->getVertexByTag(*it);
+ if(gv) {
+ getCoord(gv->x(), gv->y(), gv->z(), zeronodes[k][0],
+ zeronodes[k][1], zeronodes[k][2], gv);
_infos[k++] = AttractorInfo(*it,0,0,0);
}
- else {
- GVertex *gv = GModel::current()->getVertexByTag(*it);
- if(gv) {
- getCoord(gv->x(), gv->y(), gv->z(), zeronodes[k][0],
- zeronodes[k][1], zeronodes[k][2], gv);
- _infos[k++] = AttractorInfo(*it,0,0,0);
- }
- }
}
for(std::list<int>::iterator it = edges_id.begin();
it != edges_id.end(); ++it) {
- Curve *c = FindCurve(*it);
- if(c) {
- for(int i = 1; i < n_nodes_by_edge -1 ; i++) {
- double u = (double)i / (n_nodes_by_edge - 1);
- Vertex V = InterpolateCurve(c, u, 0);
- getCoord(V.Pos.X, V.Pos.Y, V.Pos.Z, zeronodes[k][0],
- zeronodes[k][1], zeronodes[k][2]);
+ GEdge *e = GModel::current()->getEdgeByTag(*it);
+ if(e) {
+ for(int i = 1; i < n_nodes_by_edge - 1; i++) {
+ double u = (double)i / (n_nodes_by_edge - 1);
+ Range<double> b = e->parBounds(0);
+ double t = b.low() + u * (b.high() - b.low());
+ GPoint gp = e->point(t);
+ getCoord(gp.x(), gp.y(), gp.z(), zeronodes[k][0],
+ zeronodes[k][1], zeronodes[k][2], e);
_infos[k++] = AttractorInfo(*it,1,u,0);
}
}
- else {
- GEdge *e = GModel::current()->getEdgeByTag(*it);
- if(e) {
- for(int i = 1; i < n_nodes_by_edge - 1; i++) {
- double u = (double)i / (n_nodes_by_edge - 1);
- Range<double> b = e->parBounds(0);
- double t = b.low() + u * (b.high() - b.low());
- GPoint gp = e->point(t);
- getCoord(gp.x(), gp.y(), gp.z(), zeronodes[k][0],
- zeronodes[k][1], zeronodes[k][2], e);
- _infos[k++] = AttractorInfo(*it,1,u,0);
- }
- }
- }
}
// This can lead to weird results as we generate attractors over
// the whole parametric plane (we should really use a mesh,
@@ -1745,49 +1721,37 @@ class AttractorField : public Field
int count = 0;
for(std::list<int>::iterator it = faces_id.begin();
it != faces_id.end(); ++it) {
- Surface *s = FindSurface(*it);
- if(s) {
- for(int i = 0; i < n_nodes_by_edge; i++) {
- for(int j = 0; j < n_nodes_by_edge; j++) {
- double u = (double)i / (n_nodes_by_edge - 1);
- double v = (double)j / (n_nodes_by_edge - 1);
- Vertex V = InterpolateSurface(s, u, v, 0, 0);
- getCoord(V.Pos.X, V.Pos.Y, V.Pos.Z, zeronodes[k][0],
- zeronodes[k][1], zeronodes[k][2]);
- _infos[k++] = AttractorInfo(*it,2,u,v);
- }
- }
- }
- else {
- GFace *f = GModel::current()->getFaceByTag(*it);
- if(f) {
- if (points.size()){
- for(int j = offset[count]; j < offset[count+1];j++) {
- zeronodes[k][0] = points[j].x();
- zeronodes[k][1] = points[j].y();
- zeronodes[k][2] = points[j].z();
- _infos[k++] = AttractorInfo(*it,2,uvpoints[j].x(),uvpoints[j].y());
- }
- count++;
+ GFace *f = GModel::current()->getFaceByTag(*it);
+ if(f) {
+ if (points.size()){
+ for(int j = offset[count]; j < offset[count+1];j++) {
+ zeronodes[k][0] = points[j].x();
+ zeronodes[k][1] = points[j].y();
+ zeronodes[k][2] = points[j].z();
+ _infos[k++] = AttractorInfo(*it,2,uvpoints[j].x(),uvpoints[j].y());
}
- else{
- for(int i = 0; i < n_nodes_by_edge; i++) {
- for(int j = 0; j < n_nodes_by_edge; j++) {
- double u = (double)i / (n_nodes_by_edge - 1);
- double v = (double)j / (n_nodes_by_edge - 1);
- Range<double> b1 = ge->parBounds(0);
- Range<double> b2 = ge->parBounds(1);
- double t1 = b1.low() + u * (b1.high() - b1.low());
- double t2 = b2.low() + v * (b2.high() - b2.low());
- GPoint gp = f->point(t1, t2);
- getCoord(gp.x(), gp.y(), gp.z(), zeronodes[k][0],
- zeronodes[k][1], zeronodes[k][2], f);
- _infos[k++] = AttractorInfo(*it,2,u,v);
- }
+ count++;
+ }
+ else{
+ for(int i = 0; i < n_nodes_by_edge; i++) {
+ for(int j = 0; j < n_nodes_by_edge; j++) {
+ double u = (double)i / (n_nodes_by_edge - 1);
+ double v = (double)j / (n_nodes_by_edge - 1);
+ Range<double> b1 = f->parBounds(0);
+ Range<double> b2 = f->parBounds(1);
+ double t1 = b1.low() + u * (b1.high() - b1.low());
+ double t2 = b2.low() + v * (b2.high() - b2.low());
+ GPoint gp = f->point(t1, t2);
+ getCoord(gp.x(), gp.y(), gp.z(), zeronodes[k][0],
+ zeronodes[k][1], zeronodes[k][2], f);
+ _infos[k++] = AttractorInfo(*it,2,u,v);
}
- }
- }
- }
+ }
+ }
+ }
+ else {
+ printf("face %d not yet created\n",*it);
+ }
}
kdtree = new ANNkd_tree(zeronodes, totpoints, 3);
update_needed = false;
@@ -1991,24 +1955,37 @@ void BoundaryLayerField::operator() (double x, double y, double z,
current_distance = 1.e22;
current_closest = 0;
- SMetric3 v (1./MAX_LC);
std::vector<SMetric3> hop;
+ SMetric3 v (1./(CTX::instance()->mesh.lcMax*CTX::instance()->mesh.lcMax));
+ hop.push_back(v);
for (std::list<AttractorField*>::iterator it = _att_fields.begin();
it != _att_fields.end(); ++it){
double cdist = (*(*it)) (x, y, z);
+ AttractorInfo ainfo= (*it)->getAttractorInfo().first;
SPoint3 CLOSEST= (*it)->getAttractorInfo().second;
- SMetric3 localMetric;
- if (iIntersect){
- (*this)(*it, cdist,x, y, z, localMetric, ge);
- hop.push_back(localMetric);
+ bool doNotConsider = false;
+ if (ge->dim () == ainfo.dim && ge->tag() == ainfo.ent){
+ // doNotConsider = true;
}
- if (cdist < current_distance){
- if (!iIntersect)(*this)(*it, cdist,x, y, z, localMetric, ge);
- current_distance = cdist;
- current_closest = *it;
- v = localMetric;
- _closest_point = CLOSEST;
+ else if (ge->dim () == 1 && ainfo.dim == 2){
+ GFace *gf = ge->model()->getFaceByTag(ainfo.ent);
+ // if (gf->containsEdge(ge->tag()))doNotConsider = true;
+ }
+
+ if (!doNotConsider) {
+ SMetric3 localMetric;
+ if (iIntersect){
+ (*this)(*it, cdist,x, y, z, localMetric, ge);
+ hop.push_back(localMetric);
+ }
+ if (cdist < current_distance){
+ if (!iIntersect)(*this)(*it, cdist,x, y, z, localMetric, ge);
+ current_distance = cdist;
+ current_closest = *it;
+ v = localMetric;
+ _closest_point = CLOSEST;
+ }
}
}
if (iIntersect)
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index f9f18e67617452c3d972612f30d9d56770e91df3..e1e33cdfe6bb16e700cc0bace0d37b0ba8bd712a 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -577,7 +577,18 @@ static void Mesh3D(GModel *m)
// quality reasons)
std::vector<std::vector<GRegion*> > connected;
FindConnectedRegions(delaunay, connected);
+
+ /*
for(unsigned int i = 0; i < connected.size(); i++){
+ for(unsigned j=0;j<connected[i].size();j++){
+ GRegion *gr = connected[i][j];
+ std::list<GFace*> f = gr->faces();
+ for (std::list<GFace*>::iterator it = f.begin();
+ it != f.end() ; ++it) quadsToTriangles (*it,1000000);
+ }
+ }
+ */
+ for(unsigned int i = 0; i < connected.size(); i++){
MeshDelaunayVolume(connected[i]);
//Additional code for hex mesh begin
@@ -598,7 +609,6 @@ static void Mesh3D(GModel *m)
post.execute(0);
}
}
-
}
double t2 = Cpu();
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 329c80fa33760431a97a7b1462c2784c2663fbe4..8ef3ca968eb61eb675c4524515a8e1ba5b9b82f1 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -372,7 +372,9 @@ void meshGEdge::operator() (GEdge *ge)
N = ge->meshAttributes.nbPointsTransfinite;
}
else{
- if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG || blf ){ // CTX::instance()->mesh.algo2d == ALGO_2D_PACK_PRLGRMS
+ if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG
+ // || CTX::instance()->mesh.algo2d == ALGO_2D_PACK_PRLGRMS
+ || blf){
a = Integration(ge, t_begin, t_end, F_Lc_aniso, Points,
CTX::instance()->mesh.lcIntegrationPrecision);
}
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index a877511ec008769c9fe000d1b43db1fd7c42c51f..4f8e82092c1ec5204e66719241b3d5de7596b8e6 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -305,8 +305,9 @@ static void remeshUnrecoveredEdges(std::map<MVertex*, BDS_Point*> &recoverMapInv
static bool algoDelaunay2D(GFace *gf)
{
- if(!noSeam(gf))
- return false;
+ // FIXME
+ // if(!noSeam(gf))
+ // return false;
if(gf->getMeshingAlgo() == ALGO_2D_DELAUNAY ||
gf->getMeshingAlgo() == ALGO_2D_BAMG ||
@@ -1762,6 +1763,36 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
// }
}
+ /// This is a structure that we need only for periodic cases
+ /// We will duplicate the vertices (MVertex) that are on seams
+ ///
+ /*
+ {
+ std::map<MVertex*, BDS_Point*> invertMap;
+ std::map<BDS_Point*, MVertex*>::iterator it = recoverMap.begin();
+ std::map<BDS_Point*, MVertex*>::iterator it = recoverMap.begin();
+ while(it != recoverMap.end()){
+ // we have twice vertex MVertex with 2 different coordinates
+ std::map<MVertex*, BDS_Point*>::iterator invIt = invertMap.find((*it)->second);
+ if (invIt != invertMap.end()){
+ BDS_Point* bds1 = (*invIt)->second;
+ BDS_Point* bds2 = (*it)->first;
+ MVertex * mv1 = (*it)->second;
+ // create a new "fake" vertex that will be destroyed afterwards
+ double t;
+ mv1->getParameter(0,t);
+ MVertex * mv2 = new MEdgeVertex (mv1->x(),mv1->y(),mv1->z(),mv1->onWhat(), t);
+ (*it)->second = mv2;
+ std::map<BDS_Point*, MVertex*>::iterator itTemp = it;
+ ++it;
+ recoverMap.erase(itTemp);
+ }
+ else ++it;
+ invertMap[it->second] = it->first;
+ }
+ }
+ */
+
// fill the small gmsh structures
{
std::set<BDS_Point*, PointLessThan>::iterator itp = m->points.begin();
@@ -1777,6 +1808,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
}
}
+ std::map<MTriangle*, BDS_Face*> invert_map;
{
std::list<BDS_Face*>::iterator itt = m->triangles.begin();
while (itt != m->triangles.end()){
@@ -1791,8 +1823,12 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
// when a singular point is present, degenerated triangles
// may be created, for example on a sphere that contains one
// pole
- if(v1 != v2 && v1 != v3 && v2 != v3)
+ if(v1 != v2 && v1 != v3 && v2 != v3){
+ // we are in the periodic case. if we aim at
+ // using delaunay mesh generation in thoses cases,
+ // we should double some of the vertices
gf->triangles.push_back(new MTriangle(v1, v2, v3));
+ }
}
else{
MVertex *v4 = recoverMap[n[3]];
@@ -1849,7 +1885,7 @@ void deMeshGFace::operator() (GFace *gf)
}
// for debugging, change value from -1 to -100;
-int debugSurface = -1; //-1;
+int debugSurface = -100; //-1;
void meshGFace::operator() (GFace *gf, bool print)
{
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index b603d895b8989f37333ba84c91968f7ef4a2bdf5..1f066c9b47fa04f46cec22813923d45fd99c2d39 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -1534,8 +1534,8 @@ void bowyerWatsonParallelograms(GFace *gf)
packed[i]->getParameter(1,newPoint[1]);
delete packed[i];
double metric[3];
- buildMetric(gf, newPoint, metrics[i], metric);
- //buildMetric(gf, newPoint, metric);
+ // buildMetric(gf, newPoint, metrics[i], metric);
+ buildMetric(gf, newPoint, metric);
bool success = insertAPoint( gf, AllTris.begin(), newPoint, metric, DATA , AllTris, 0, oneNewTriangle, &oneNewTriangle);
if (!success) oneNewTriangle = 0;
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 7d119ce2dd09ad3c4a5b5ac20032c5008896bef3..b5bb52269d883b728df408c488716d1e8f059506 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -2182,7 +2182,7 @@ void laplaceSmoothing(GFace *gf, int niter, bool infinity_norm)
int untangleInvalidQuads(GFace *gf, int niter)
{
- // return;
+ // return 0;
int N = 0;
#if defined(HAVE_BFGS)
v2t_cont adj;
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 14916f6f1a4d44d78be28212128aac2a3eab1d38..3c9fc1b91e9630a152c9060b346f421d8ed1277c 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -1031,6 +1031,7 @@ void meshNormalsPointOutOfTheRegion(GRegion *gr)
void meshGRegion::operator() (GRegion *gr)
{
+
gr->model()->setCurrentMeshEntity(gr);
if(gr->geomType() == GEntity::DiscreteVolume) return;
diff --git a/Mesh/surfaceFiller.cpp b/Mesh/surfaceFiller.cpp
index 371ebc628c20ac4b97be44c2e8ada54da9c33b16..f20fcf562d49bd9d83fdfe96b9e1ab18503ff028 100644
--- a/Mesh/surfaceFiller.cpp
+++ b/Mesh/surfaceFiller.cpp
@@ -18,8 +18,9 @@
#include "BackgroundMesh.h"
#include "intersectCurveSurface.h"
-static const double FACTOR = .71;
+static const double FACTOR = .81;
static const int NUMDIR = 3;
+//static const double DIRS [NUMDIR] = {0.0};
static const double DIRS [NUMDIR] = {0.0, M_PI/20.,-M_PI/20.};
/// a rectangle in the tangent plane is transformed
@@ -35,6 +36,7 @@ struct surfacePointWithExclusionRegion {
SPoint2 _p[4][NUMDIR];
SPoint2 _q[4];
SMetric3 _meshMetric;
+ double _distanceSummed;
/*
+ p3
p4 |
@@ -44,12 +46,23 @@ struct surfacePointWithExclusionRegion {
*/
- surfacePointWithExclusionRegion (MVertex *v, SPoint2 p[4][NUMDIR], SMetric3 & meshMetric){
- _meshMetric = meshMetric;
+ surfacePointWithExclusionRegion (MVertex *v, SPoint2 p[4][NUMDIR], SMetric3 & meshMetric, surfacePointWithExclusionRegion *father = 0){
_v = v;
+ _meshMetric = meshMetric;
_center = (p[0][0]+p[1][0]+p[2][0]+p[3][0])*.25;
for (int i=0;i<4;i++)_q[i] = _center + (p[i][0]+p[(i+1)%4][0]-_center*2)*FACTOR;
for (int i=0;i<4;i++)for (int j=0;j<NUMDIR;j++)_p[i][j] = p[i][j];
+
+ if (!father){
+ fullMatrix<double> V(3,3);
+ fullVector<double> S(3);
+ meshMetric.eig(V,S);
+ double l = std::max(std::max(S(0),S(1)),S(2));
+ _distanceSummed = sqrt(1/(l*l));
+ }
+ else {
+ _distanceSummed = father->_distanceSummed + distance (father->_v,_v);
+ }
}
bool inExclusionZone (const SPoint2 &p){
double mat[2][2];
@@ -88,6 +101,19 @@ struct my_wrapper {
my_wrapper (SPoint2 sp) : _tooclose (false), _p(sp) {}
};
+class compareSurfacePointWithExclusionRegionPtr
+{
+ public:
+ inline bool operator () (const surfacePointWithExclusionRegion *a, const surfacePointWithExclusionRegion *b) const
+ {
+ if(a->_distanceSummed > b->_distanceSummed) return false;
+ if(a->_distanceSummed < b->_distanceSummed) return true;
+ return a<b;
+ }
+};
+
+
+
bool rtree_callback(surfacePointWithExclusionRegion *neighbour,void* point){
my_wrapper *w = static_cast<my_wrapper*>(point);
@@ -294,7 +320,8 @@ void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, std::vec
}
// put boundary vertices in a fifo queue
- std::queue<surfacePointWithExclusionRegion*> fifo;
+ // std::queue<surfacePointWithExclusionRegion*> fifo;
+ std::set<surfacePointWithExclusionRegion*, compareSurfacePointWithExclusionRegionPtr> fifo;
std::vector<surfacePointWithExclusionRegion*> vertices;
// put the RTREE
RTree<surfacePointWithExclusionRegion*,double,2,double> rtree;
@@ -305,7 +332,8 @@ void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, std::vec
compute4neighbors (gf, *it, false, newp, metricField);
surfacePointWithExclusionRegion *sp =
new surfacePointWithExclusionRegion (*it, newp, metricField);
- fifo.push(sp);
+ // fifo.push(sp);
+ fifo.insert(sp);
vertices.push_back(sp);
double _min[2],_max[2];
sp->minmax(_min,_max);
@@ -319,8 +347,10 @@ void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, std::vec
while(!fifo.empty()){
//surfacePointWithExclusionRegion & parent = fifo.top();
- surfacePointWithExclusionRegion * parent = fifo.front();
- fifo.pop();
+ // surfacePointWithExclusionRegion * parent = fifo.front();
+ surfacePointWithExclusionRegion * parent = *fifo.begin();
+ // fifo.pop();
+ fifo.erase(fifo.begin());
for (int dir=0;dir<NUMDIR;dir++){
// printf("dir = %d\n",dir);
int countOK = 0;
@@ -335,8 +365,9 @@ void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, std::vec
// printf(" %g %g %g %g\n",parent._center.x(),parent._center.y(),gp.u(),gp.v());
compute4neighbors (gf, v, false, newp, metricField);
surfacePointWithExclusionRegion *sp =
- new surfacePointWithExclusionRegion (v, newp, metricField);
- fifo.push(sp);
+ new surfacePointWithExclusionRegion (v, newp, metricField, parent);
+ // fifo.push(sp);
+ fifo.insert(sp);
vertices.push_back(sp);
double _min[2],_max[2];
sp->minmax(_min,_max);
diff --git a/Post/PView.cpp b/Post/PView.cpp
index a1920dab2be4b55b1a57d4ec660648cad9b77651..30da150af0f7b0b78db1a4d166db7c0548a1c6b5 100644
--- a/Post/PView.cpp
+++ b/Post/PView.cpp
@@ -121,6 +121,8 @@ PView::PView(const std::string &name, const std::string &type,
t = PViewDataGModel::ElementData;
else if(type == "ElementNodeData")
t = PViewDataGModel::ElementNodeData;
+ else if(type == "Beam")
+ t = PViewDataGModel::BeamData;
else{
Msg::Error("Unknown type of view to create '%s'", type.c_str());
return;
diff --git a/Post/PViewDataGModel.h b/Post/PViewDataGModel.h
index 08d21619196f409c4eba49dd8c2fd10db86cd3d1..869078e0d08457f58f426c7b56f4d80a2a23cb6a 100644
--- a/Post/PViewDataGModel.h
+++ b/Post/PViewDataGModel.h
@@ -154,7 +154,8 @@ class PViewDataGModel : public PViewData {
NodeData = 1,
ElementData = 2,
ElementNodeData = 3,
- GaussPointData = 4
+ GaussPointData = 4,
+ BeamData = 5
};
private:
// the data, indexed by time step
diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index 6cd11d435202d9ce040027f8c351a41bf81b9e09..36cfcd598098b469ca967acfc717416907c83c0f 100644
--- a/Solver/CMakeLists.txt
+++ b/Solver/CMakeLists.txt
@@ -16,6 +16,7 @@ elasticitySolver.cpp
multiscaleLaplace.cpp
functionSpace.cpp
filters.cpp
+ frameSolver.cpp
sparsityPattern.cpp
STensor43.cpp
STensor33.cpp
diff --git a/Solver/frameSolver.cpp b/Solver/frameSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..294a9ec1085505fd54b1dbb0f6091e0e21a9a262
--- /dev/null
+++ b/Solver/frameSolver.cpp
@@ -0,0 +1,288 @@
+#include "GmshConfig.h"
+#include "GModel.h"
+#include "GVertex.h"
+#include "GEdge.h"
+#include "frameSolver.h"
+#include "linearSystemCSR.h"
+#include "linearSystemPETSc.h"
+#include "linearSystemGMM.h"
+#include "linearSystemFull.h"
+#if defined(HAVE_POST)
+#include "PView.h"
+#include "PViewData.h"
+#endif
+
+frameSolver2d::frameSolver2d (GModel *gm)
+ : pAssembler(0),_myModel(gm)
+{
+}
+
+
+void frameSolver2d::addFixations (const std::vector<int> & dirs, const std::vector<int> &modelVertices, double value)
+{
+ for (unsigned int j=0;j<modelVertices.size();j++){
+ GVertex *gv = _myModel->getVertexByTag(modelVertices[j]);
+ if (gv){
+ for (unsigned int i=0;i<dirs.size();i++){
+ _fixations.push_back(gmshFixation(gv,dirs[i],value));
+ }
+ }
+ }
+}
+
+void frameSolver2d::addNodalForces (const std::vector<int> &modelVertices, const std::vector<double> & force)
+{
+ for (unsigned int j=0;j<modelVertices.size();j++){
+ GVertex *gv = _myModel->getVertexByTag(modelVertices[j]);
+ if (gv){
+ _nodalForces.push_back(std::make_pair(gv,force));
+ }
+ }
+}
+
+void frameSolver2d::addBeamsOrBars (const std::vector<int> &modelEdges,
+ double E, double I, double A, int r[2]){
+ int r_middle[2] = {1,1} , r_left[2] = {r[0],1} , r_right[2] = {0,r[1]};
+ // printf("adding %d beams\n",modelEdges.size());
+ for (unsigned int i=0;i<modelEdges.size();i++){
+ GEdge *ge = _myModel->getEdgeByTag(modelEdges[i]);
+ if (ge){
+ // printf("model edge %d found\n",ge->tag());
+ for (unsigned int j=0; j < ge->lines.size(); ++j){
+ MLine *l = ge->lines[j];
+ if (j == 0 && j == ge->lines.size()-1) _beams.push_back(gmshBeam2d ( l , E , I , A , r));
+ else if (j == 0) _beams.push_back(gmshBeam2d ( l , E , I , A , r_left));
+ else if (j == ge->lines.size()-1) _beams.push_back(gmshBeam2d ( l , E , I , A , r_right));
+ else _beams.push_back(gmshBeam2d ( l , E , I , A , r_middle));
+ }
+ }
+ }
+}
+
+void frameSolver2d::addBeams (const std::vector<int> &modelEdges,
+ double E, double I, double A)
+{
+ int r[2] = {1,1};
+ addBeamsOrBars(modelEdges,E,I,A,r);
+}
+
+void frameSolver2d::addBars (const std::vector<int> &modelEdges,
+ double E, double I, double A)
+{
+ int r[2] = {0,0};
+ addBeamsOrBars(modelEdges,E,I,A,r);
+}
+
+
+
+// solve any time a parameter is modified
+/*
+
+ A vertex is connected to beams. The question
+ is how many bars are rotuled
+
+ We define 2 dofs per node
+
+ */
+void frameSolver2d::createDofs ()
+{
+ // printf("coucou %d fixations\n",_fixations.size());
+ for (unsigned int i = 0; i<_fixations.size(); ++i){
+ const gmshFixation &f = _fixations[i];
+ // printf("f._vertex(%d) = %p %d %g\n",i,f._vertex,f._direction,f._value);
+ MVertex *v = f._vertex->mesh_vertices[0];
+ Dof DOF (v->getNum(),f._direction);
+ pAssembler->fixDof(DOF, f._value);
+ }
+
+ // printf("coucou2\n");
+ computeRotationTags ();
+ // printf("coucou3\n");
+ for (unsigned int i=0;i<_beams.size(); i++){
+ // printf("beam[%d] rot %d %d\n",i,_beams[i]._rotationTags[0],_beams[i]._rotationTags[1]);
+ for (unsigned int j=0;j<2;j++){
+ MVertex *v = _beams[i]._element->getVertex(j);
+ Dof theta (v->getNum(),
+ Dof::createTypeWithTwoInts (2,_beams[i]._rotationTags[j]));
+ pAssembler->numberDof(theta);
+ Dof U (v->getNum(),0);
+ pAssembler->numberDof(U);
+ Dof V (v->getNum(),1);
+ pAssembler->numberDof(V);
+ }
+ }
+ // printf("%d dofs\n",pAssembler->sizeOfR());
+}
+
+void frameSolver2d::computeStiffnessMatrix (int iBeam,
+ fullMatrix<double> &K)
+{
+ const gmshBeam2d &b = _beams[iBeam];
+ const double BS = b._E * b._I / (b._L*b._L*b._L);
+ const double TS = b._E * b._A / b._L;
+ const MVertex *v1 = b._element->getVertex(0);
+ const MVertex *v2 = b._element->getVertex(1);
+ const double alpha = atan2 (v2->y() - v1->y(),
+ v2->x() - v1->x());
+ const double C = cos(alpha);
+ const double S = sin(alpha);
+
+ printf("beam %d %g %g %g\n",iBeam,alpha,C,S);
+
+ fullMatrix<double> R (6,6);
+ R(0,0) = R(1,1) = R (3,3) = R(4,4) = C;
+ R(0,1) = R (3,4) = S;
+ R(1,0) = R (4,3) = -S;
+ R(2,2) = R(5,5) = 1.0;
+
+ fullMatrix<double> k (6,6);
+
+ // tensile stiffness
+ k(0,0) = k (3,3) = TS;
+ k(0,3) = k (3,0) = -TS;
+ // bending stiffness
+ k(1,1) = k(4,4) = 12 * BS;
+ k(2,2) = k(5,5) = 4. * BS * b._L * b._L;
+ k(1,2) = k(2,1) = k(1,5) = k(5,1) = 6 * BS * b._L;
+ k(4,2) = k(2,4) = k(4,5) = k(5,4) = -6 * BS * b._L;
+ k(4,1) = k(1,4) = -12 * BS ;
+ k(5,2) = k(2,5) = -2 * BS * b._L * b._L;
+
+ fullMatrix<double> Rt (R), temp(6,6);
+ Rt.transposeInPlace();
+ Rt.mult(k,temp);temp.mult(R,K);
+}
+
+void frameSolver2d::solve ()
+{
+#if defined(HAVE_TAUCS)
+ linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
+#elif defined(HAVE_PETSC)
+ linearSystemPETSc<double> *lsys = new linearSystemPETSc<double>;
+#else
+ linearSystemGmm<double> *lsys = new linearSystemGmm<double>;
+ lsys->setNoisy(2);
+#endif
+
+ if (pAssembler)delete pAssembler;
+ pAssembler = new dofManager<double>(lsys);
+
+ // fix dofs and create free ones
+ createDofs();
+
+ // force vector
+ std::vector<std::pair<GVertex*,std::vector<double> > >::iterator it =
+ _nodalForces.begin();
+ for (; it != _nodalForces.end(); ++it){
+ MVertex *v = it->first->mesh_vertices[0];
+ const std::vector<double> & F = it->second;
+ Dof DOFX (v->getNum(),0);
+ Dof DOFY (v->getNum(),1);
+ pAssembler->assemble(DOFX, F[0]);
+ pAssembler->assemble(DOFY, F[1]);
+ }
+
+ // stifness matrix
+ for (unsigned int i=0;i<_beams.size(); i++){
+ fullMatrix<double> K(6,6);
+ computeStiffnessMatrix (i,K);
+ _beams[i]._stiffness = K;
+ MVertex *v0 = _beams[i]._element->getVertex(0);
+ MVertex *v1 = _beams[i]._element->getVertex(1);
+ Dof theta0 (v0->getNum(),
+ Dof::createTypeWithTwoInts (2,_beams[i]._rotationTags[0]));
+ Dof theta1 (v1->getNum(),
+ Dof::createTypeWithTwoInts (2,_beams[i]._rotationTags[1]));
+ Dof U0 (v0->getNum(),0);
+ Dof U1 (v1->getNum(),0);
+ Dof V0 (v0->getNum(),1);
+ Dof V1 (v1->getNum(),1);
+ Dof DOFS[6] = {U0,V0,theta0,U1,V1,theta1};
+ for (int j=0;j<6;j++){
+ for (int k=0;k<6;k++){
+ pAssembler->assemble(DOFS[j],DOFS[k],K(j,k));
+ }
+ }
+ }
+ lsys->systemSolve();
+
+ // save the solution
+ for (unsigned int i=0;i<_beams.size(); i++){
+ MVertex *v0 = _beams[i]._element->getVertex(0);
+ MVertex *v1 = _beams[i]._element->getVertex(1);
+ Dof theta0 (v0->getNum(),
+ Dof::createTypeWithTwoInts (2,_beams[i]._rotationTags[0]));
+ Dof theta1 (v1->getNum(),
+ Dof::createTypeWithTwoInts (2,_beams[i]._rotationTags[1]));
+ Dof U0 (v0->getNum(),0);
+ Dof U1 (v1->getNum(),0);
+ Dof V0 (v0->getNum(),1);
+ Dof V1 (v1->getNum(),1);
+ Dof DOFS[6] = {U0,V0,theta0,U1,V1,theta1};
+ for (int j=0;j<6;j++){
+ pAssembler->getDofValue(DOFS[j], _beams[i]._displacement[j]);
+ }
+ }
+ delete lsys;
+ delete pAssembler;
+}
+
+void frameSolver2d::exportFrameData(const char * DISPL, const char * M) {
+ {
+ std::map<int, std::vector<double> > data;
+ for (unsigned int i=0;i<_beams.size(); i++){
+ std::vector<double> tmp;
+ // tmp.push_back(_beams[i]._E);
+ // tmp.push_back(_beams[i]._I);
+ // tmp.push_back(_beams[i]._A);
+ for (int j=0;j<6;j++){
+ tmp.push_back(_beams[i]._displacement[j]);
+ }
+ data[_beams[i]._element->getNum()] = tmp;
+ }
+ PView *pv = new PView ("displacements", "Beam", _myModel, data, 0.0, 6);
+ pv->getData()->writeMSH(DISPL);
+ delete pv;
+ }
+ {
+ std::map<int, std::vector<double> > data;
+ for (unsigned int i=0;i<_beams.size(); i++){
+ std::vector<double> tmp;
+ fullVector<double> d (_beams[i]._displacement,6), F(6);
+ _beams[i]._stiffness.mult(d,F);
+ tmp.push_back(-F(2));
+ tmp.push_back(F(5));
+ data[_beams[i]._element->getNum()] = tmp;
+ }
+ PView *pv = new PView ("Momentum", "ElementNodeData", _myModel, data, 0.0, 1);
+ pv->getData()->writeMSH(M);
+ delete pv;
+ }
+
+}
+
+void frameSolver2d::computeRotationTags ()
+{
+ std::multimap<MVertex*,gmshBeam2d*> v2b;
+ for (unsigned int i=0;i<_beams.size(); i++){
+ v2b.insert(std::make_pair(_beams[i]._element->getVertex(0),&_beams[i]));
+ v2b.insert(std::make_pair(_beams[i]._element->getVertex(1),&_beams[i]));
+ }
+
+ std::multimap<MVertex*,gmshBeam2d*>::iterator s_it;
+ for (std::multimap<MVertex*,gmshBeam2d*>::iterator it = v2b.begin();
+ it != v2b.end(); it = s_it){
+ MVertex *theKey = it->first;
+
+ std::pair<std::multimap<MVertex*,gmshBeam2d*>::iterator,
+ std::multimap<MVertex*,gmshBeam2d*>::iterator>
+ keyRange = v2b.equal_range(theKey);
+ int countRotules = 0;
+ for (s_it = keyRange.first; s_it != keyRange.second; ++s_it){
+ gmshBeam2d *b = s_it->second;
+ if (!b->isRigid(theKey)){
+ b->setRotationTag(theKey,++countRotules);
+ }
+ }
+ }
+}
diff --git a/Solver/frameSolver.h b/Solver/frameSolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..9bf9ea753bfb771fca81bc66f8449e224ffcca0c
--- /dev/null
+++ b/Solver/frameSolver.h
@@ -0,0 +1,77 @@
+// Gmsh - Copyright (C) 1997-2013 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to the public mailing list <gmsh@geuz.org>.
+
+#ifndef _FRAME_SOLVER_H_
+#define _FRAME_SOLVER_H_
+
+class GVertex;
+
+#include <map>
+#include <string>
+#include "PView.h"
+#include "dofManager.h"
+#include "fullMatrix.h"
+#include "MLine.h"
+
+struct gmshBeam2d
+{
+ MLine *_element;
+ double _I, _A , _E, _L;
+ bool _rigidNodes [2];
+ double _forceVector[6];
+ double _displacement[6];
+ int _rotationTags [2];
+ fullMatrix<double> _stiffness;
+ gmshBeam2d (MLine *l, double E , double I, double A, int r[2])
+ : _element(l), _A(A), _E(E), _I(I){
+ _L = distance(_element->getVertex(0),_element->getVertex(1));
+ _rotationTags[0] = _rotationTags[1] = 0;
+ _rigidNodes[0] = r[0];_rigidNodes[1] = r[1];
+ for (int i=0;i<6;i++)_displacement[i] = 0.;
+ }
+ inline bool isRigid(MVertex*v)const{
+ return _element->getVertex(0) == v ? _rigidNodes[0] : _rigidNodes[1];
+ }
+ inline void setRotationTag(MVertex * v, int tag){
+ _element->getVertex(0) == v ? _rotationTags[0] = tag : _rotationTags[1] = tag;
+ }
+};
+
+struct gmshFixation
+{
+ GVertex * _vertex;
+ int _direction;
+ double _value;
+ gmshFixation (GVertex *v, int d, double val) : _vertex(v),_direction(d), _value(val) {}
+};
+
+// a solver that computes equilibriums of frames...
+class frameSolver2d
+{
+ dofManager<double> *pAssembler;
+ std::vector<gmshBeam2d> _beams;
+ std::vector<std::pair<GVertex*,std::vector<double> > > _nodalForces;
+ std::vector<gmshFixation > _fixations;
+ GModel *_myModel;
+ void computeStiffnessMatrix (int iBeam, fullMatrix<double> &K);
+ void createDofs () ;
+ void computeRotationTags () ;
+ void addBeamsOrBars (const std::vector<int> &modelEdges,
+ double E, double I, double A, int r[2]);
+
+
+ public:
+ frameSolver2d (GModel* myModel);
+ void addBeams (const std::vector<int> &modelEdges,
+ double E, double I, double A);
+ void addBars (const std::vector<int> &modelEdges,
+ double E, double I, double A);
+ void addNodalForces (const std::vector<int> &modelVertices, const std::vector<double> & force);
+ void addFixations (const std::vector<int> & dirs, const std::vector<int> &modelVertices, double value);
+ void exportFrameData(const char * displ, const char* M);
+ void solve () ;
+};
+
+#endif
diff --git a/benchmarks/2d/conge.geo b/benchmarks/2d/conge.geo
index 1b8384cccd3e57a9866d9babb8713e98e6b3be2e..39841f15cba0f75e96eb5196ef566bfedfac0625 100644
--- a/benchmarks/2d/conge.geo
+++ b/benchmarks/2d/conge.geo
@@ -78,5 +78,5 @@ Physical Line(25) = {12,13,14,15,16,17,18,19,20,10};
Physical Surface(26) = {24,22};
Physical Line(27) = {10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 14, 13, 12, 11};
Physical Line(28) = {17, 16, 20, 19, 18, 15};
-//Recombine Surface {24, 22};
+Recombine Surface {24, 22};
Mesh.RecombinationAlgorithm=1;
\ No newline at end of file
diff --git a/wrappers/gmshpy/gmshSolver.i b/wrappers/gmshpy/gmshSolver.i
index e4e815a3ccb3d4c2e64dbdadf2198b57185216c7..61e1a2da181ac8e045fccd0d7ab4dd4ab7560207 100644
--- a/wrappers/gmshpy/gmshSolver.i
+++ b/wrappers/gmshpy/gmshSolver.i
@@ -8,6 +8,7 @@
#if defined(HAVE_SOLVER)
#include "dofManager.h"
#include "elasticitySolver.h"
+ #include "frameSolver.h"
#include "linearSystem.h"
#include "linearSystemCSR.h"
#include "linearSystemFull.h"
@@ -20,6 +21,7 @@
%include "dofManager.h"
%template(dofManagerDouble) dofManager<double>;
%include "elasticitySolver.h"
+%include "frameSolver.h"
%include "linearSystem.h"
%template(linearSystemDouble) linearSystem<double>;
%template(linearSystemFullMatrixDouble) linearSystem<fullMatrix<double> >;