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

clean up + fix bug in recurElement (gbricteux)

parent e3eea890
Branches
Tags
No related merge requests found
...@@ -36,7 +36,8 @@ PView* elasticitySolver::buildDisplacementView (const std::string &postFileName ...@@ -36,7 +36,8 @@ PView* elasticitySolver::buildDisplacementView (const std::string &postFileName
std::vector<double> d; std::vector<double> d;
d.push_back(0); d.push_back(0); d.push_back(0); d.push_back(0); d.push_back(0); d.push_back(0);
for (int i = 0; i < elasticFields.size(); ++i){ for (int i = 0; i < elasticFields.size(); ++i){
const double E = (elasticFields[i]._enrichment) ? (*elasticFields[i]._enrichment)((*it)->x(),(*it)->y(),(*it)->z()) : 1.0; const double E = (elasticFields[i]._enrichment) ?
(*elasticFields[i]._enrichment)((*it)->x(), (*it)->y(), (*it)->z()) : 1.0;
// printf("%g %d \n",pAssembler->getDofValue(*it,0,elasticFields[i]._tag),elasticFields[i]._tag); // printf("%g %d \n",pAssembler->getDofValue(*it,0,elasticFields[i]._tag),elasticFields[i]._tag);
d[0] += pAssembler->getDofValue(*it, 0, elasticFields[i]._tag) * E; d[0] += pAssembler->getDofValue(*it, 0, elasticFields[i]._tag) * E;
d[1] += pAssembler->getDofValue(*it, 1, elasticFields[i]._tag) * E; d[1] += pAssembler->getDofValue(*it, 1, elasticFields[i]._tag) * E;
...@@ -176,6 +177,10 @@ void elasticitySolver::readInputFile(const std::string &fn) ...@@ -176,6 +177,10 @@ void elasticitySolver::readInputFile(const std::string &fn)
if(fscanf(f, "%s", name) != 1) return; if(fscanf(f, "%s", name) != 1) return;
setMesh(name); setMesh(name);
} }
else {
Msg::Error("Invalid input : %s", what);
return;
}
} }
fclose(f); fclose(f);
} }
......
...@@ -105,7 +105,6 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const ...@@ -105,7 +105,6 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
} }
printf("\n"); printf("\n");
printf("\n"); printf("\n");
} }
void elasticityTerm::elementVector(SElement *se, fullVector<double> &m) const void elasticityTerm::elementVector(SElement *se, fullVector<double> &m) const
......
This diff is collapsed.
...@@ -37,13 +37,13 @@ void recurCut(RecurElement *re, int maxlevel, int level) ...@@ -37,13 +37,13 @@ void recurCut(RecurElement *re, int maxlevel, int level)
RecurElement *re0 = new RecurElement(&DI_Triangle(p1, p12, p13)); RecurElement *re0 = new RecurElement(&DI_Triangle(p1, p12, p13));
recurCut(re0, maxlevel, level); recurCut(re0, maxlevel, level);
re->sub[0] = re0; re0->super = re; re->sub[0] = re0; re0->super = re;
RecurElement *re1 = new RecurElement(&DI_Triangle(p2, p12, p23)); RecurElement *re1 = new RecurElement(&DI_Triangle(p2, p23, p12));
recurCut(re1, maxlevel, level); recurCut(re1, maxlevel, level);
re->sub[1] = re1; re1->super = re; re->sub[1] = re1; re1->super = re;
RecurElement *re2 = new RecurElement(&DI_Triangle(p3, p13, p23)); RecurElement *re2 = new RecurElement(&DI_Triangle(p3, p13, p23));
recurCut(re2, maxlevel, level); recurCut(re2, maxlevel, level);
re->sub[2] = re2; re2->super = re; re->sub[2] = re2; re2->super = re;
RecurElement *re3 = new RecurElement(&DI_Triangle(p12, p13, p23)); RecurElement *re3 = new RecurElement(&DI_Triangle(p12, p23, p13));
recurCut(re3, maxlevel, level); recurCut(re3, maxlevel, level);
re->sub[3] = re3; re3->super = re; re->sub[3] = re3; re3->super = re;
} }
...@@ -87,22 +87,22 @@ void recurCut(RecurElement *re, int maxlevel, int level) ...@@ -87,22 +87,22 @@ void recurCut(RecurElement *re, int maxlevel, int level)
RecurElement *re0 = new RecurElement(&DI_Tetra(p1, p12, p13, p14)); RecurElement *re0 = new RecurElement(&DI_Tetra(p1, p12, p13, p14));
recurCut(re0, maxlevel, level); recurCut(re0, maxlevel, level);
re->sub[0] = re0; re0->super = re; re->sub[0] = re0; re0->super = re;
RecurElement *re1 = new RecurElement(&DI_Tetra(p2, p12, p23, p24)); RecurElement *re1 = new RecurElement(&DI_Tetra(p2, p23, p12, p24));
recurCut(re1, maxlevel, level); recurCut(re1, maxlevel, level);
re->sub[1] = re1; re1->super = re; re->sub[1] = re1; re1->super = re;
RecurElement *re2 = new RecurElement(&DI_Tetra(p3, p13, p23, p34)); RecurElement *re2 = new RecurElement(&DI_Tetra(p3, p13, p23, p34));
recurCut(re2, maxlevel, level); recurCut(re2, maxlevel, level);
re->sub[2] = re2; re2->super = re; re->sub[2] = re2; re2->super = re;
RecurElement *re3 = new RecurElement(&DI_Tetra(p4, p14, p24, p34)); RecurElement *re3 = new RecurElement(&DI_Tetra(p4, p14, p34, p24));
recurCut(re3, maxlevel, level); recurCut(re3, maxlevel, level);
re->sub[3] = re3; re3->super = re; re->sub[3] = re3; re3->super = re;
RecurElement *re4 = new RecurElement(&DI_Tetra(p12, p14, p24, p34)); RecurElement *re4 = new RecurElement(&DI_Tetra(p12, p14, p24, p34));
recurCut(re4, maxlevel, level); recurCut(re4, maxlevel, level);
re->sub[4] = re4; re4->super = re; re->sub[4] = re4; re4->super = re;
RecurElement *re5 = new RecurElement(&DI_Tetra(p12, p23, p24, p34)); RecurElement *re5 = new RecurElement(&DI_Tetra(p12, p23, p34, p24));
recurCut(re5, maxlevel, level); recurCut(re5, maxlevel, level);
re->sub[5] = re5; re5->super = re; re->sub[5] = re5; re5->super = re;
RecurElement *re6 = new RecurElement(&DI_Tetra(p12, p13, p23, p34)); RecurElement *re6 = new RecurElement(&DI_Tetra(p12, p13, p34, p23));
recurCut(re6, maxlevel, level); recurCut(re6, maxlevel, level);
re->sub[6] = re6; re6->super = re; re->sub[6] = re6; re6->super = re;
RecurElement *re7 = new RecurElement(&DI_Tetra(p12, p13, p14, p34)); RecurElement *re7 = new RecurElement(&DI_Tetra(p12, p13, p14, p34));
...@@ -140,25 +140,25 @@ void recurCut(RecurElement *re, int maxlevel, int level) ...@@ -140,25 +140,25 @@ void recurCut(RecurElement *re, int maxlevel, int level)
RecurElement *re0 = new RecurElement(&DI_Hexa(p1, p12, p1234, p14, p15, p1256, p12345678, p1458)); RecurElement *re0 = new RecurElement(&DI_Hexa(p1, p12, p1234, p14, p15, p1256, p12345678, p1458));
recurCut(re0, maxlevel, level); recurCut(re0, maxlevel, level);
re->sub[0] = re0; re0->super = re; re->sub[0] = re0; re0->super = re;
RecurElement *re1 = new RecurElement(&DI_Hexa(p2, p23, p1234, p12, p26, p2367, p12345678, p1256)); RecurElement *re1 = new RecurElement(&DI_Hexa(p12, p2, p23, p1234, p1256, p26, p2367, p12345678));
recurCut(re1, maxlevel, level); recurCut(re1, maxlevel, level);
re->sub[1] = re1; re1->super = re; re->sub[1] = re1; re1->super = re;
RecurElement *re2 = new RecurElement(&DI_Hexa(p3, p34, p1234, p23, p37, p3478, p12345678, p2367)); RecurElement *re2 = new RecurElement(&DI_Hexa(p1234, p23, p3, p34, p12345678, p2367, p37, p3478));
recurCut(re2, maxlevel, level); recurCut(re2, maxlevel, level);
re->sub[2] = re2; re2->super = re; re->sub[2] = re2; re2->super = re;
RecurElement *re3 = new RecurElement(&DI_Hexa(p4, p14, p1234, p34, p48, p1458, p12345678, p3478)); RecurElement *re3 = new RecurElement(&DI_Hexa(p14, p1234, p34, p4, p1458, p12345678, p3478, p48));
recurCut(re3, maxlevel, level); recurCut(re3, maxlevel, level);
re->sub[3] = re3; re3->super = re; re->sub[3] = re3; re3->super = re;
RecurElement *re4 = new RecurElement(&DI_Hexa(p5, p58, p5678, p56, p15, p1458, p12345678, p1256)); RecurElement *re4 = new RecurElement(&DI_Hexa(p15, p1256, p12345678, p1458, p5, p56, p5678, p58));
recurCut(re4, maxlevel, level); recurCut(re4, maxlevel, level);
re->sub[4] = re4; re4->super = re; re->sub[4] = re4; re4->super = re;
RecurElement *re5 = new RecurElement(&DI_Hexa(p6, p56, p5678, p67, p26, p1256, p12345678, p2367)); RecurElement *re5 = new RecurElement(&DI_Hexa(p1256, p26, p2367, p12345678, p56, p6, p67, p5678));
recurCut(re5, maxlevel, level); recurCut(re5, maxlevel, level);
re->sub[5] = re5; re5->super = re; re->sub[5] = re5; re5->super = re;
RecurElement *re6 = new RecurElement(&DI_Hexa(p7, p67, p5678, p78, p37, p2367, p12345678, p3478)); RecurElement *re6 = new RecurElement(&DI_Hexa(p12345678, p2367, p37, p3478, p5678, p67, p7, p78));
recurCut(re6, maxlevel, level); recurCut(re6, maxlevel, level);
re->sub[6] = re6; re6->super = re; re->sub[6] = re6; re6->super = re;
RecurElement *re7 = new RecurElement(&DI_Hexa(p8, p78, p5678, p58, p48, p3478, p12345678, p1458)); RecurElement *re7 = new RecurElement(&DI_Hexa(p1458, p12345678, p3478, p48, p58, p5678, p78, p8));
recurCut(re7, maxlevel, level); recurCut(re7, maxlevel, level);
re->sub[7] = re7; re7->super = re; re->sub[7] = re7; re7->super = re;
} }
...@@ -186,7 +186,7 @@ bool signChange (RecurElement *re, const DI_Element *e, const std::vector<const ...@@ -186,7 +186,7 @@ bool signChange (RecurElement *re, const DI_Element *e, const std::vector<const
else { else {
for(unsigned int p = 0; p < cp.size(); p++) for(unsigned int p = 0; p < cp.size(); p++)
cp[p].chooseLs(Lsi); cp[p].chooseLs(Lsi);
re->el->chooseLs(Lsi); if (re->super) re->el->chooseLs(Lsi);
} }
} }
re->el->clearLs(); re->el->clearLs();
......
...@@ -7,7 +7,7 @@ OCC = /Users/remacle/SOURCES/opencascade/lib/libTKSTEP.a /Users/remacle/SOURCES/ ...@@ -7,7 +7,7 @@ OCC = /Users/remacle/SOURCES/opencascade/lib/libTKSTEP.a /Users/remacle/SOURCES/
mainElasticity: mainElasticity.cpp mainElasticity: mainElasticity.cpp
g++ ${FLAGS} -o mainElasticity ${INC} mainElasticity.cpp\ g++ ${FLAGS} -o mainElasticity ${INC} mainElasticity.cpp\
${GMSH} -llapack -lblas -lm -L/Users/remacle/SOURCES/gmsh/contrib/taucs/lib -ltaucs ${GMSH} -llapack -lblas -lm -lmetis -L/home/gaetan/develop/gmsh/contrib/taucs/lib/linux -ltaucs
mainCartesian: mainCartesian.cpp mainCartesian: mainCartesian.cpp
g++ ${FLAGS} -o mainCartesian ${INC} mainCartesian.cpp\ g++ ${FLAGS} -o mainCartesian ${INC} mainCartesian.cpp\
...@@ -31,7 +31,7 @@ mainPost: mainPost.cpp ...@@ -31,7 +31,7 @@ mainPost: mainPost.cpp
mainLevelset: mainLevelset.cpp mainLevelset: mainLevelset.cpp
g++ ${FLAGS} -o mainLevelset ${INC} mainLevelset.cpp\ g++ ${FLAGS} -o mainLevelset ${INC} mainLevelset.cpp\
${GMSH} -llapack -lblas ${GMSH} -llapack -lblas -lm -lmetis -L/home/gaetan/develop/gmsh/contrib/taucs/lib/linux -ltaucs
clean: clean:
rm -rf mainSimple mainAntTweakBar mainGlut mainPost mainLevelset\ rm -rf mainSimple mainAntTweakBar mainGlut mainPost mainLevelset\
......
...@@ -12,6 +12,7 @@ int main (int argc, char* argv[]){ ...@@ -12,6 +12,7 @@ int main (int argc, char* argv[]){
// globals are still present in Gmsh // globals are still present in Gmsh
GmshInitialize(argc, argv); GmshInitialize(argc, argv);
GmshSetOption("General","Terminal",1.);
// instanciate a solver // instanciate a solver
elasticitySolver mySolver (1000); elasticitySolver mySolver (1000);
...@@ -20,7 +21,7 @@ int main (int argc, char* argv[]){ ...@@ -20,7 +21,7 @@ int main (int argc, char* argv[]){
mySolver.readInputFile(argv[1]); mySolver.readInputFile(argv[1]);
// solve the problem // solve the problem
mySolver.solve(); mySolver.solve(); printf("problem solved\n");
PView *pv = mySolver.buildDisplacementView("displacement"); PView *pv = mySolver.buildDisplacementView("displacement");
pv->getData()->writeMSH("disp.msh", false); pv->getData()->writeMSH("disp.msh", false);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment