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

*** empty log message ***

parent e30d06c7
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Geo/GFaceCompound.cpp
+ 15
− 9
View file @ fee07df2
......@@ -11,6 +11,7 @@
#include "Octree.h"
#include "gmshAssembler.h"
#include "gmshLaplace.h"
#include "gmshConvexCombination.h"
#include "gmshLinearSystemGmm.h"
#include "gmshLinearSystemFull.h"
#include "FunctionSpace.h"
......@@ -50,12 +51,14 @@ public:
class gmshDistanceBasedDiffusivity : public gmshFunction<double>
{
private:
mutable int comp;
MElement *_current;
mutable std::map<MVertex*, SPoint3> _coordinates;
public:
gmshDistanceBasedDiffusivity (std::map<MVertex*,SPoint3> &coordinates)
: _current (0),_coordinates(coordinates){}
: _current (0),_coordinates(coordinates), comp(0){}
void setCurrent (MElement *current){ _current = current; }
void setCompound(int compound){ comp = compound; }
virtual double operator () (double x, double y, double z) const
{
double xyz[3] = {x, y, z}, uvw[3];
......@@ -66,8 +69,11 @@ public:
value[i] = p[2];
}
double val = _current->interpolate(value, uvw[0], uvw[1], uvw[2]);
//return exp(5*val);
//return 1./(exp(x)+1.e-4);//exp(5*val);
return 1.0;
//printf("compiound =%d \n", comp);
//if (comp == 8) return 1.0;
//else return 1.e-15;
}
};
......@@ -98,13 +104,9 @@ void GFaceCompound::getBoundingEdges()
for (std::list<GFace*>::iterator it = _compound.begin(); it != _compound.end(); ++it){
printf("set compound %d for face %d \n", tag(),(*it)->tag());
(*it)->setCompound(this);
}
printf("***** In GFaceCompound: size U0=%d, v0=%d\n ", _U0.size(), _V0.size());
//in case the bounding edges are explicitely given
if (_U0.size()){
std::list<GEdge*> :: const_iterator it = _U0.begin();
......@@ -146,7 +148,7 @@ void GFaceCompound::getBoundingEdges()
std::set<GEdge*>::iterator itf = _unique.begin();
for ( ; itf != _unique.end(); ++itf){
l_edges.push_back(*itf);
printf("for edge %d, add face %d \n", (*itf)->tag(), this->tag());
//printf("for edge %d, add face %d \n", (*itf)->tag(), this->tag());
(*itf)->addFace(this);
}
......@@ -491,12 +493,15 @@ void GFaceCompound::parametrize(iterationStep step) const
}
}
else{
//gmshConvexCombinationTerm laplace(model(), &ONE, 1);
gmshLaplaceTerm laplace(model(), &ONE, 1);
//gmshLaplaceTerm laplace(model(), &diffusivity, 1);
it = _compound.begin();
for ( ; it != _compound.end() ; ++it){
for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
MTriangle *t = (*it)->triangles[i];
diffusivity.setCurrent(t);
//diffusivity.setCompound((*it)->tag());
//diffusivity.setCurrent(t);
laplace.addToMatrix(myAssembler, t);
}
}
......@@ -917,7 +922,8 @@ void GFaceCompound::printStuff() const
t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
K1, K2, K3);
it0->second.z(),it1->second.z(),it2->second.z());
//K1, K2, K3);
fprintf(uvx,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
it0->second.x(), it0->second.y(), 0.0,
......
Geo/GModel.cpp
+ 5
− 7
View file @ fee07df2
......@@ -916,7 +916,7 @@ int GModel::removeDuplicateMeshVertices(double tolerance)
void GModel::createTopologyFromMesh()
{
printf("***** In createTopologyFromMesh: \n");
Msg::Debug("***** In createTopologyFromMesh:");
std::vector<GEntity*> entities;
getEntities(entities);
......@@ -944,10 +944,10 @@ void GModel::createTopologyFromMesh()
}
}
printf("vertices size =%d \n", Dvertices.size());
printf("edges size =%d \n", Dedges.size());
printf("faces size =%d \n", Dfaces.size());
printf("regions size =%d \n", Dregions.size());
Msg::Debug("vertices size =%d", Dvertices.size());
Msg::Debug("edges size =%d ", Dedges.size());
Msg::Debug("faces size =%d ", Dfaces.size());
Msg::Debug("regions size =%d", Dregions.size());
//For each discreteRegion, create Topology
//---------------------------------------
......@@ -1094,7 +1094,6 @@ void GModel::createTopologyFromMesh()
vE = v1;
i=-1;
}
if (it == myEdges.end()) break;
}
printf("end Edges \n");
......@@ -1132,7 +1131,6 @@ void GModel::createTopologyFromMesh()
if (std::find(all_vertices.begin(), all_vertices.end(), v1) == all_vertices.end()) all_vertices.push_back(v1);
}
e->mesh_vertices.insert(e->mesh_vertices.begin(), all_vertices.begin(), all_vertices.end());
printf("all vertice size =%d\n", all_vertices.size());
for (std::vector<int>::iterator itFace = tagFaces.begin(); itFace != tagFaces.end(); itFace++) {
GFace *dFace = getFaceByTag(abs(*itFace));
......
Geo/GRegionCompound.cpp
+ 3
− 3
View file @ fee07df2
......@@ -15,7 +15,7 @@
GRegionCompound::GRegionCompound(GModel *m, int tag, std::vector<GRegion*> &compound)
: GRegion(m, tag), _compound(compound)
{
printf("********* In GRegion compound size =%d \n", _compound.size());
//printf("********* In GRegion compound size =%d \n", _compound.size());
getBoundingFaces();
}
......@@ -29,7 +29,7 @@ void GRegionCompound::getBoundingFaces(){
std::multiset<GFace*> _touched;
std::vector<GRegion*>::iterator it = _compound.begin();
for ( ; it != _compound.end(); ++it){
printf("face %d \n", (*it)->tag());
//printf("face %d \n", (*it)->tag());
std::list<GFace*> ed = (*it)->faces();
std::list<GFace*> :: iterator ite = ed.begin();
for ( ; ite != ed.end(); ++ite){
......@@ -49,7 +49,7 @@ void GRegionCompound::getBoundingFaces(){
std::set<GFace*>::iterator itf = _unique.begin();
for ( ; itf != _unique.end(); ++itf){
l_faces.push_back(*itf);
printf("for face %d, add region %d \n", (*itf)->tag(), this->tag());
//printf("for face %d, add region %d \n", (*itf)->tag(), this->tag());
(*itf)->addRegion(this);
}
......
Geo/Makefile
+ 2
− 0
View file @ fee07df2
......@@ -137,6 +137,8 @@ GFaceCompound${OBJEXT}: GFaceCompound.cpp ../Common/GmshConfig.h GFaceCompound.h
../Geo/GEntity.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h \
../Geo/MElement.h ../Numeric/gmshFunction.h ../Common/Gmsh.h \
../Common/GmshMessage.h ../Numeric/GmshMatrix.h \
../Numeric/gmshConvexCombination.h ../Numeric/gmshTermOfFormulation.h \
../Numeric/gmshFunction.h ../Numeric/GmshMatrix.h \
../Numeric/gmshLinearSystemGmm.h ../Numeric/gmshLinearSystem.h \
../Numeric/gmshLinearSystemFull.h ../Numeric/gmshLinearSystem.h \
../Numeric/GmshMatrix.h
......
Geo/discreteEdge.cpp
+ 0
− 11
View file @ fee07df2
......@@ -338,15 +338,6 @@ void discreteEdge::parametrize()
}
// for (int i = 0; i < mesh_vertices.size(); i++){
// double t1;
// mesh_vertices[i]->getParameter(0,t1);
// printf("** AFTER v1=%d t1=%g\n", mesh_vertices[i]->getNum(),t1 );
// }
computeNormals();
}
......@@ -387,8 +378,6 @@ void discreteEdge::computeNormals () const
//printf("normal =%g %g %g \n", itn->second.x(), itn->second.y(), itn->second.z());
}
printf("********* normal size = %d \n", _normals.size());
}
void discreteEdge::getLocalParameter(const double &t, int &iLine,
......
Numeric/Makefile
+ 19
− 0
View file @ fee07df2
......@@ -14,6 +14,7 @@ CFLAGS = ${OPTIM} ${FLAGS} ${INC} ${SYSINCLUDE}
SRC = Numeric.cpp\
GmshMatrix.cpp\
gmshLaplace.cpp\
gmshConvexCombination.cpp\
gmshHelmholtz.cpp\
gmshElasticity.cpp\
gmshProjection.cpp\
......@@ -79,6 +80,24 @@ gmshLaplace${OBJEXT}: gmshLaplace.cpp gmshLaplace.h gmshTermOfFormulation.h \
../Geo/MEdge.h ../Geo/SVector3.h ../Numeric/FunctionSpace.h \
../Numeric/GmshMatrix.h ../Numeric/Gauss.h ../Common/Gmsh.h \
../Common/GmshMessage.h Numeric.h
gmshConvexCombination${OBJEXT}: gmshConvexCombination.cpp \
gmshConvexCombination.h gmshTermOfFormulation.h GmshMatrix.h \
../Common/GmshConfig.h ../Common/GmshMessage.h gmshFunction.h \
gmshAssembler.h gmshLinearSystem.h ../Geo/GModel.h ../Geo/GVertex.h \
../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/SOrientedBoundingBox.h \
../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/Pair.h \
../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h \
../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/GEdgeLoop.h \
../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Common/GmshDefines.h \
../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
../Geo/MEdge.h ../Geo/SVector3.h ../Numeric/FunctionSpace.h \
../Numeric/GmshMatrix.h ../Numeric/Gauss.h ../Common/Gmsh.h \
../Common/GmshMessage.h Numeric.h
gmshHelmholtz${OBJEXT}: gmshHelmholtz.cpp gmshHelmholtz.h \
gmshTermOfFormulation.h GmshMatrix.h ../Common/GmshConfig.h \
../Common/GmshMessage.h gmshFunction.h gmshAssembler.h \
......
Numeric/gmshConvexCombination.cpp 0 → 100644
+ 28
− 0
View file @ fee07df2
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
#include "gmshConvexCombination.h"
#include "Numeric.h"
void gmshConvexCombinationTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
{
m.set_all(0.);
int nbNodes = e->getNumVertices();
double x = 0.3*(e->getVertex(0)->x()+e->getVertex(1)->x()+e->getVertex(2)->x());
double y = 0.3*(e->getVertex(0)->y()+e->getVertex(1)->y()+e->getVertex(2)->y());
double z = 0.3*(e->getVertex(0)->z()+e->getVertex(1)->z()+e->getVertex(2)->z());
const double _diff = 1.0;
//const double _diff = (*_diffusivity)(x, y, z);
for (int j = 0; j < nbNodes; j++){
for (int k = 0; k < nbNodes; k++){
m(j,k) = -1.*_diff;
}
m(j,j) = (nbNodes-1)*_diff;
}
}
Numeric/gmshConvexCombination.h 0 → 100644
+ 36
− 0
View file @ fee07df2
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
#ifndef _GMSH_CONVEX_H_
#define _GMSH_CONVEX_H_
#include "gmshTermOfFormulation.h"
#include "gmshFunction.h"
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
#include "GmshMatrix.h"
class gmshConvexCombinationTerm : public gmshNodalFemTerm<double> {
protected:
const gmshFunction<double> *_diffusivity;
const int _iField ;
public:
virtual int sizeOfR(MElement *e) const { return e->getNumVertices(); }
virtual int sizeOfC(MElement *e) const { return e->getNumVertices(); }
void getLocalDofR(MElement *e, int iRow, MVertex **vR, int *iCompR, int *iFieldR) const
{
*vR = e->getVertex(iRow);
*iCompR = 0;
*iFieldR = _iField;
}
public:
gmshConvexCombinationTerm(GModel *gm, gmshFunction<double> *diffusivity, int iField = 0) :
gmshNodalFemTerm<double>(gm), _diffusivity(diffusivity), _iField(iField){}
virtual void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
};
#endif
Numeric/gmshLaplace.cpp
+ 27
− 7
View file @ fee07df2
......@@ -38,14 +38,14 @@ void gmshLaplaceTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
invjac[2][2] * grads[j][2];
}
double pi = 3.14;
double H_x=1.0;
double H_y=1.0;
double H_z=1.0;
double K_x=1.0;
double K_y=1.0;
double K_z=1.0;
for (int j = 0; j < nbNodes; j++){
for (int k = 0; k <= j; k++){
m(j, k) += (H_x*Grads[j][0] * Grads[k][0] +
H_y*Grads[j][1] * Grads[k][1] +
H_z*Grads[j][2] * Grads[k][2]) * weight * detJ * _diff;
m(j, k) += (K_x*Grads[j][0] * Grads[k][0] +
K_y*Grads[j][1] * Grads[k][1] +
K_z*Grads[j][2] * Grads[k][2]) * weight * detJ * _diff;
}
}
}
......@@ -53,6 +53,26 @@ void gmshLaplaceTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
for (int j = 0; j < nbNodes; j++)
for (int k = 0; k < j; k++)
m(k, j) = m(j, k);
//check for positive scheme
// //printf("ELEM\n");
for (int j = 0; j < nbNodes; j++){
double sum = m(j,0)+m(j,1)+m(j,2);
if (sum < 0.0){
printf("Sum LINE gmshLaplace Term NEG <0 %g %g %g\n", m(j,0),m(j,1),m(j,2));
// for (int k = 0; k < nbNodes; k++) m(j,k) = -1.;
// m(j,j) = (nbNodes-1);
}
// else{
// printf("Sum POS\n");
// }
// //for (int k = 0; k < nbNodes; k++) {
// //printf("m(%d,%d)=%g ", j,k, m(j,k));
// //}
// //printf("\n");
}
}
void gmshLaplaceTerm2DParametric::elementMatrix(MElement *e, gmshMatrix<double> &m) const
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment