|
LR-splines
0.5
|
#include <cstdio>#include <cmath>#include <cstdlib>#include <iostream>#include <fstream>#include <cstring>#include "LRSpline/LRSplineSurface.h"#include "LRSpline/LRSplineVolume.h"#include "LRSpline/Element.h"#include "LRSpline/Basisfunction.h"Go to the source code of this file.
Functions | |
| int | main (int argc, char **argv) |
| int main | ( | int | argc, |
| char ** | argv | ||
| ) |
Definition at line 17 of file Diagonal.cpp.
References LR::LRSpline::generateIDs(), LR::LRSpline::getAllBasisfunctions(), LR::LRSpline::getAllElements(), LR::LRSplineSurface::getDiagonalBasisfunctions(), LR::LRSplineSurface::getDiagonalElements(), LR::LRSplineVolume::getDiagonalElements(), LR::LRSplineSurface::isLinearIndepByMappingMatrix(), LR_FULLSPAN, LR_MINSPAN, LR_STRUCTURED_MESH, LR::LRSpline::nBasisFunctions(), LR::LRSpline::nElements(), LR::LRSplineSurface::nMeshlines(), LR::LRSplineVolume::nMeshRectangles(), LR::LRSplineSurface::refineBasisFunction(), LR::LRSplineSurface::refineElement(), LR::LRSplineVolume::refineElement(), LR::LRSpline::setRefMultiplicity(), LR::LRSpline::setRefStrat(), LR::LRSplineSurface::writePostscriptElements(), LR::LRSplineSurface::writePostscriptFunctionSpace(), and LR::LRSplineSurface::writePostscriptMesh().
{
// set default parameter values
int p = 3;
int N = 6;
int m = 1;
int scheme = 0;
bool dumpfile = false;
bool vol = false;
string parameters(" parameters: \n" \
" -p <n> polynomial DEGREE (order-1) of the basis\n" \
" -n <n> number of iterations\n" \
" -m <n> knot multiplicity\n" \
" -dumpfile create .eps and .m files of the meshes\n" \
" -vol enforce a volumetric test case\n"\
" -scheme <n> refinement scheme (0=FULLSPAN, 1=MINSPAN, 2=STRUCT)\n" \
" -help display (this) help information\n");
// read input
for(int i=1; i<argc; i++) {
if(strcmp(argv[i], "-p") == 0)
p = atoi(argv[++i]);
else if(strcmp(argv[i], "-scheme") == 0)
scheme = atoi(argv[++i]);
else if(strcmp(argv[i], "-n") == 0)
N = atoi(argv[++i]);
else if(strcmp(argv[i], "-m") == 0)
m = atoi(argv[++i]);
else if(strcmp(argv[i], "-dumpfile") == 0)
dumpfile = true;
else if(strcmp(argv[i], "-vol") == 0)
vol = true;
else if(strcmp(argv[i], "-help") == 0) {
cout << "usage: " << argv[0] << endl << parameters;
exit(0);
} else {
cout << "usage: " << argv[0] << endl << parameters;
exit(1);
}
}
// do some error testing on input
if(m > p )
m = p;
if(scheme > 2)
scheme = 2;
// clear previous result from file
if(dumpfile) {
ofstream clear_out;
clear_out.open("dof.m");
clear_out.close();
clear_out.open("elements.m");
clear_out.close();
}
// make two identical splines
LRSplineSurface *lrs;
LRSplineVolume *lrv;
if(vol) {
lrv = new LRSplineVolume(p+1, p+1, p+1, p+1, p+1, p+1);
} else {
lrs = new LRSplineSurface(p+1, p+1, p+1, p+1);
}
// setup refinement parameters
enum refinementStrategy strat;
if(scheme == 0)
strat = LR_FULLSPAN;
else if(scheme == 1)
strat = LR_MINSPAN;
else if(scheme == 2)
strat = LR_STRUCTURED_MESH;
if(vol) {
lrv->setRefMultiplicity(m);
lrv->setRefStrat(strat);
} else {
lrs->setRefMultiplicity(m);
lrs->setRefStrat(strat);
}
// for all iterations
for(int n=0; n<N; n++) {
vector<int> indices;
if(scheme < 2) {
if(vol) {
lrv->generateIDs();
lrv->getDiagonalElements(indices);
lrv->refineElement(indices);
} else {
lrs->generateIDs();
lrs->getDiagonalElements(indices);
lrs->refineElement(indices);
}
} else if(scheme == 2) {
if(vol) {
// lrv->generateIDs();
// lrv->getDiagonalBasisfunctions(indices);
// lrv->refineBasisFunction(indices);
} else {
lrs->generateIDs();
lrs->getDiagonalBasisfunctions(indices);
lrs->refineBasisFunction(indices);
}
}
// draw result files (with next refinement-step-diagonal shaded)
if(dumpfile && !vol) {
vector<int> diagElms, diagFuncs;
lrs->getDiagonalElements(diagElms);
lrs->getDiagonalBasisfunctions(diagFuncs);
char filename[128];
sprintf(filename, "index_%02d.eps", n+1);
ofstream out;
out.open(filename);
if(scheme < 2)
lrs->writePostscriptMesh(out, true, &diagElms);
else
lrs->writePostscriptMesh(out);
out.close();
sprintf(filename, "parameter_%02d.eps", n+1);
out.open(filename);
if(scheme < 2)
lrs->writePostscriptElements(out, 2,2, true, &diagElms);
else
lrs->writePostscriptElements(out, 2,2, true);
out.close();
sprintf(filename, "func_%02d.eps", n+1);
out.open(filename);
if(scheme == 2)
lrs->writePostscriptFunctionSpace(out, &diagFuncs);
else
lrs->writePostscriptFunctionSpace(out);
out.close();
out.open("dof.m", ios::app);
out << lrs->nBasisFunctions() << endl;
out.close();
out.open("elements.m", ios::app);
out << lrs->nElements() << endl;
out.close();
}
}
vector<int> overloadedBasis;
vector<int> overloadedElements;
vector<int> multipleOverloadedElements;
// harvest some statistics and display these results
int nBasis = 0;
int nElements = 0;
int nMeshlines = 0;
if(vol) lrv->generateIDs();
else lrs->generateIDs();
double avgBasisToElement = 0;
double avgBasisToLine = 0;
int maxBasisToElement = -1;
int minBasisToElement = 9999999;
int nOverloadedElms = 0;
int nOverloadedBasis = 0;
if(vol) {
nBasis = lrv->nBasisFunctions();
nElements = lrv->nElements();
nMeshlines = lrv->nMeshRectangles();
} else {
nBasis = lrs->nBasisFunctions();
nElements = lrs->nElements();
nMeshlines = lrs->nMeshlines();
}
HashSet<Basisfunction*> basis = (vol) ? lrv->getAllBasisfunctions() : lrs->getAllBasisfunctions();
vector<Element*> elements = (vol) ? lrv->getAllElements() : lrs->getAllElements();
for(Basisfunction* b : basis) {
int nE = b->nSupportedElements();
maxBasisToElement = (maxBasisToElement > nE) ? maxBasisToElement : nE;
minBasisToElement = (minBasisToElement < nE) ? minBasisToElement : nE;
avgBasisToElement += nE;
if(b->isOverloaded()) {
nOverloadedBasis++;
overloadedBasis.push_back(b->getId());
}
}
avgBasisToElement /= nBasis;
avgBasisToLine /= nBasis;
double avgElementToBasis = 0;
double avgSquareElementToBasis = 0;
double hashCodePercentage = ((double)basis.uniqueHashCodes())/nBasis;
int maxElementToBasis = -1;
int minElementToBasis = 9999999;
for(Element *e : elements) {
int nB = e->nBasisFunctions();
maxElementToBasis = (maxElementToBasis > nB) ? maxElementToBasis : nB;
minElementToBasis = (minElementToBasis < nB) ? minElementToBasis : nB;
avgElementToBasis += nB;
avgSquareElementToBasis += nB*nB;
if(e->isOverloaded()) {
nOverloadedElms++;
overloadedElements.push_back(e->getId());
}
}
avgElementToBasis /= nElements;
avgSquareElementToBasis /= nElements;
cout << "Some statistics: " << endl;
cout << "-------------------------------------------------------------" << endl;
cout << "Number of basisfunctions: " << nBasis << endl;
cout << "Number of elements : " << nElements << endl;
cout << "Number of meshlines : " << nMeshlines << endl;
cout << "-------------------------------------------------------------" << endl;
cout << "Min number of Basisfuntion -> Element: " << minBasisToElement << endl;
cout << "Max number of Basisfuntion -> Element: " << maxBasisToElement << endl;
cout << "Avg number of Basisfuntion -> Element: " << avgBasisToElement << endl;
cout << endl;
cout << "Min number of Element -> Basisfunction: " << minElementToBasis << endl;
cout << "Max number of Element -> Basisfunction: " << maxElementToBasis << endl;
cout << "Avg number of Element -> Basisfunction: " << avgElementToBasis << endl;
cout << "Avg square number of Element -> Basisfunction: " << avgSquareElementToBasis
<< " (" << sqrt(avgSquareElementToBasis) << ")" << endl;
cout << endl;
cout << "Number of overloaded Basisfunctions : " << nOverloadedBasis << endl;
cout << "Number of overloaded Elements : " << nOverloadedElms << endl;
cout << endl;
cout << "Number of unique hashcodes : " << basis.uniqueHashCodes() ;
cout << " (" << hashCodePercentage*100 << " %)" << endl;
cout << "-------------------------------------------------------------" << endl;
if(nBasis < 1300 && !vol) {
cout << "Is linearly independent : " << ((lrs->isLinearIndepByMappingMatrix(false) )? "True":"False") << endl;
cout << "-------------------------------------------------------------" << endl;
}
if(dumpfile) {
cout << endl;
cout << "Written";
if(!vol) {
ofstream meshfile;
meshfile.open("mesh.eps");
lrs->writePostscriptMesh(meshfile);
meshfile.close();
cout << " mesh to mesh.eps and";
}
ofstream lrfile;
lrfile.open("diagonal.lr");
if(vol) lrfile << *lrv << endl;
else lrfile << *lrs << endl;
lrfile.close();
cout << " diagonal.lr\n";
}
}
1.7.6.1