wzsxgrid.cxx
来自「有限元学习研究用源代码(老外的),供科研人员参考」· CXX 代码 · 共 1,773 行 · 第 1/3 页
CXX
1,773 行
sx->priv->FEMFactors=FEM2Edge;
sx->priv->sxNext=sxNextFace2;
sx->priv->faceiter2= new DListIter<EDG2>(sx->grid->mesh->castToMESH2()->lines());
break;
case 3:
sx->nlines=3;
sx->priv->FVMFactors=FVM3Triangle;
sx->priv->FEMFactors=FEM3Triangle;
sx->priv->sxNext=sxNextFace3;
sx->priv->faceiter3= new DListIter<TR3>(sx->grid->mesh->castToMESH3()->triangles());
break;
}
sx->dim=g->spacedim-1;
sx->codim=1;
break;
case SX_EDGES:
sx->nnodes=2;
sx->nlines=1;
sx->linenodes=0;
switch(sx->spacedim)
{
case 1:
sx->priv->FVMFactors=0;
sx->priv->FEMFactors=0;
sx->priv->sxNext=sxNextLine1;
sx->priv->lineiter1= new DListIter<EDG1>(sx->grid->mesh->castToMESH1()->lines());
break;
case 2:
sx->priv->FVMFactors=0;
sx->priv->FEMFactors=0;
sx->priv->sxNext=sxNextLine2;
sx->priv->lineiter2= new DListIter<EDG2>(sx->grid->mesh->castToMESH2()->lines());
break;
case 3:
sx->priv->FVMFactors=0;
sx->priv->FEMFactors=0;
sx->priv->sxNext=sxNextLine3;
sx->priv->lineiter3= new DListIter<EDG3>(sx->grid->mesh->castToMESH3()->lines());
break;
}
sx->dim=1;
sx->codim=g->spacedim-1;
break;
case SX_NODES:
sx->nnodes=1; /* hier koennten Nachbarknoten/kanten stehen */
sx->nlines=0;
sx->linenodes=NULL;
sx->node=node_node_ptr;
sx->dim=0;
sx->codim=g->spacedim;
switch(sx->spacedim)
{
case 1:
sx->priv->sxNext=sxNextNode1;
sx->priv->nodeiter1= new DListIter<PT1>(sx->grid->mesh->castToMESH1()->points());
break;
case 2:
sx->priv->sxNext=sxNextNode2;
sx->priv->nodeiter2= new DListIter<PT2>(sx->grid->mesh->castToMESH2()->points());
break;
case 3:
sx->priv->sxNext=sxNextNode3;
sx->priv->nodeiter3= new DListIter<PT3>(sx->grid->mesh->castToMESH3()->points());
break;
}
sx->priv->FVMFactors=FVM0Vertex;
sx->priv->FEMFactors=FEM0Vertex;
break;
default:
dprintf(err)("sorry, loop number %d not implemented\n",WhatToTraverse);
}
{
int i,j;
sx->stiffmat[0]=0;
sx->massmat[0]=0;
for (i=1,j=1;i<=sx->nnodes;i++,j+=sx->nnodes)
sx->stiffmat[i]=&sx->priv->stiffmat[j]-1,
sx->massmat[i]=&sx->priv->massmat[j]-1;
}
for(in=1;in<=sx->nnodes;in++)
sx->coord[in]=sx->priv->coord+(in-1)*sx->spacedim;
sx->priv->sxNext(sx);
}
int sxStop(Simplex sx)
{
return (sx->priv==0);
}
////////////////////////////////////////////////////
// various data
void* sxMesh(sxGrid g) { return g->mesh; }
int sxDimension(sxGrid g) { return g->spacedim; }
int sxNumber(sxGrid g, int comp)
{
switch (comp)
{
case SX_NODES: return g->node_num;
case SX_CELLS: return g->cell_num;
case SX_FACES: return g->face_num;
case SX_EDGES: return g->line_num;
case SX_ALL_FACES: return g->all_face_num;
default: return 0;
}
}
int sxNumberOfMaterials(sxGrid g, int comp)
{
switch (comp)
{
case SX_NODES: return g->node_mat_num;
case SX_CELLS: return g->cell_mat_num;
case SX_FACES: return g->face_mat_num;
case SX_EDGES: return 0;
default: return 0;
}
}
double sxXMin(sxGrid g) { return g->xmin; }
double sxXMax(sxGrid g) { return g->xmax; }
double sxYMin(sxGrid g) { return g->ymin; }
double sxYMax(sxGrid g) { return g->ymax; }
double sxZMin(sxGrid g) { return g->zmin; }
double sxZMax(sxGrid g) { return g->zmax; }
CmdPars Cmd; // command Parameters
//extern - on SGI initialization seems to work only in main prog!
////////////////////////////////////////////////////
// Create/Destroy
///////////////////////////////////////////////////
void sxgPrintData(sxGrid g) // same for creation and refinement!
{
int i;
int bad;
double *cell_mat_volume;
SimplexData sx;
int bad_cell_count_zero_vol;
int bad_cell_count_neg_line;
double bad_cell_volume;
bad_cell_count_zero_vol=0;
bad_cell_count_neg_line=0;
bad_cell_volume=0;
cell_mat_volume=(double*)calloc(g->cell_mat_num+1,sizeof(double));
SX_LOOP(g,sx,SX_CELLS)
{
bad=0;
sxFVMFactors(&sx);
cell_mat_volume[0]+=sx.volume;
cell_mat_volume[sx.material]+=sx.volume;
if (sx.volume==0.0)
{
bad_cell_count_zero_vol++;
bad=1;
}
for (i=1;i<=sx.nlines;i++)
{
if (sx.linefac[i]<0.0)
{
bad_cell_count_neg_line++;
bad=1;
break;
}
}
if (bad) bad_cell_volume+=sx.volume;
}
if (debug(inf))
printf("-----------------------------------------------------------------\n");
dprintf(inf)("cell data:\n");
dprintf(inf)("grid: %4d cells, vol=%g\n",g->cell_num,cell_mat_volume[0]);
for(i=1;i<=g->cell_mat_num;i++)
if(g->cell_mat_cell_num[i])
dprintf(inf)("material %d: %4d cells, vol=%g\n",
i,
g->cell_mat_cell_num[i],cell_mat_volume[i]);
dprintf(inf)("bad cells: zero vol: %d, neg line: %d, vol=%g\n",
bad_cell_count_zero_vol,
bad_cell_count_neg_line,
bad_cell_volume);
free(cell_mat_volume);
if (debug(inf))
printf("-----------------------------------------------------------------\n");
dprintf(inf)("face data:\n");
for(i=1;i<=g->face_mat_num;i++)
if(g->face_mat_face_num[i])
dprintf(inf)("material %d: %d faces\n",i,g->face_mat_face_num[i]);
dprintf(inf)("boundary faces: %d\n",g->face_num);
dprintf(inf)("all faces: %d\n",g->all_face_num);
if (debug(inf))
printf("-----------------------------------------------------------------\n");
dprintf(inf)("node data:\n");
dprintf(inf)("nodes: %d\n",g->node_num);
dprintf(inf)("spacedim: %d\n",g->spacedim);
dprintf(inf)("xcoord: [ %g , %g ]\n",g->xmin,g->xmax);
dprintf(inf)("ycoord: [ %g , %g ]\n",g->ymin,g->ymax);
if (g->spacedim>2)
dprintf(inf)("zcoord: [ %g , %g ]\n",g->zmin,g->zmax);
if (debug(inf))
printf("-----------------------------------------------------------------\n");
}
static int bisect=0;
void sxSetBisection(void)
{
bisect=1;
}
void number_stuff(sxGrid g)
{
MESH *m=g->mesh;
register int i;
// sorry for this!
switch (g->spacedim)
{
case 1:
{
PT1 *patch;
DListIter<PT1> iter(m->castToMESH1()->points());
i=0;
while ((patch = iter.all()))
{
i++;
patch->setNode(i);
}
g->node_num=i;
}
break;
case 2:
{
PT2 *patch;
DListIter<PT2> iter(m->castToMESH2()->points());
i=0;
while ((patch = iter.all()))
{
i++;
patch->setNode(i);
}
g->node_num=i;
}
break;
case 3:
{
PT3 *patch;
DListIter<PT3> iter(m->castToMESH3()->points());
i=0;
while ((patch =iter.all()))
{
i++;
patch->setNode(i++);
}
g->node_num=i;
}
break;
}
PATCH *patch;
m->resetElemIter();
for(i=1;i<=g->cell_mat_num;i++) g->cell_mat_cell_num[i]=0;
i=0;
while ((patch = m->elemIterAll()))
{
i++;
patch->setNode(i);
g->cell_mat_cell_num[patch->Class()]++;
}
switch(g->spacedim)
{
case 1:
{
EDG1* patch;
DListIter<EDG1> iter(m->castToMESH1()->lines());
i=0;
while ((patch = iter.all()))
{
i++;
patch->setNode(i);
}
g->line_num=i;
}
break;
case 2:
{
EDG2* patch;
DListIter<EDG2> iter(m->castToMESH2()->lines());
i=0;
while ((patch = iter.all()))
{
i++;
patch->setNode(i);
}
g->line_num=i;
}
break;
case 3:
{
EDG3* patch;
DListIter<EDG3> iter(m->castToMESH3()->lines());
i=0;
while ((patch = iter.all()))
{
i++;
patch->setNode(i);
}
g->line_num=i;
}
break;
}
for(i=1;i<=g->face_mat_num;i++) g->face_mat_face_num[i]=0;
switch(g->spacedim)
{
case 1:
{
PT1* patch;
DListIter<PT1> iter(m->castToMESH1()->points());
while ((patch = iter.all()))
if (patch->onBoundary()&&patch->Class())
g->face_mat_face_num[patch->Class()]++;
g->all_face_num=g->line_num;
}
break;
case 2:
{
EDG2* patch;
DListIter<EDG2> iter(m->castToMESH2()->lines());
while ((patch = iter.all()))
if (patch->onBoundary()&&patch->Class())
g->face_mat_face_num[patch->Class()]++;
g->all_face_num=g->line_num;
}
break;
case 3:
{
TR3* patch;
DListIter<TR3> iter(m->castToMESH3()->triangles());;
iter.reset();
int i=0;
while ((patch = iter.all()))
{
i++;
patch->setNode(i);
if (patch->onBoundary()&&patch->Class())
g->face_mat_face_num[patch->Class()]++;
}
g->all_face_num=i;
}
break;
}
g->face_num=0;
for(i=1;i<=g->face_mat_num;i++) g->face_num+=g->face_mat_face_num[i];
g->cell_num=g->mesh->noOfElements();
/*
{
void find_face_in_list(n1,n2,n3)
{
sort(n1,n2,n3);
}
}
*/
}
sxGrid sxCreate(char *meshname, void *user_data)
{
register int i;
sxGrid g;
char fname[64];
MESH *m=NULL;
int spacedim=0;
init();
Cmd.ReadCmdFile(meshname);
///////////////////////////////////////
///// MESH creation
if (Cmd.isTrue("printParameters")) cout << Cmd;
if (spacedim==0)
{
if (!Cmd.get("spaceDim",&spacedim))
{
dprintf(wrn)("Spacedim said to be 0, assuming 2!\n");
spacedim=2;
}
}
sprintf(fname,"%s.geo",meshname);
switch (spacedim)
{
case 1: m=new MESH1(fname);
break;
case 2:
if (bisect)
{
m=new MESH2BISECT(fname);
}
else
m=new MESH2(fname);
break;
case 3: m=new MESH3(fname);
break;
}
assert(m!=NULL);
g=new sxGridData;
assert(g!=NULL);
g->mesh=m;
g->spacedim=spacedim;
if (bisect)
{
g->bisect=bisect;
bisect=0;
}
///////////////////////////////////////
///// Set Simplex format data
switch (g->spacedim)
{
case 1:
g->cell_line_num=1;
g->cell_face_num=2;
g->cell_line_node_ptr[1][1]=1;
g->cell_line_node_ptr[1][2]=2;
g->cell_face_node_ptr[1][1]=1;
g->cell_face_node_ptr[1][2]=2;
g->cell_face_node_ptr[2][1]=2;
g->cell_face_node_ptr[2][2]=1;
break;
case 2:
g->cell_line_num=3;
g->cell_face_num=3;
g->cell_line_node_ptr[1][1]=2;
g->cell_line_node_ptr[1][2]=3;
g->cell_line_node_ptr[1][3]=1;
g->cell_line_node_ptr[2][1]=1;
g->cell_line_node_ptr[2][2]=3;
g->cell_line_node_ptr[2][3]=2;
g->cell_line_node_ptr[3][1]=1;
g->cell_line_node_ptr[3][2]=2;
g->cell_line_node_ptr[3][3]=3;
g->cell_face_node_ptr[1][1]=2;
g->cell_face_node_ptr[1][2]=3;
g->cell_face_node_ptr[1][3]=1;
g->cell_face_node_ptr[2][1]=1;
g->cell_face_node_ptr[2][2]=3;
g->cell_face_node_ptr[2][3]=2;
g->cell_face_node_ptr[3][1]=1;
g->cell_face_node_ptr[3][2]=2;
g->cell_face_node_ptr[3][3]=3;
break;
case 3:
g->cell_line_num=6;
g->cell_face_num=4;
/* die ersten beiden muessen mit celltype konform gehen !*/
g->cell_line_node_ptr[1][1]=3;
g->cell_line_node_ptr[1][2]=4;
g->cell_line_node_ptr[1][3]=1;
g->cell_line_node_ptr[1][4]=2;
g->cell_line_node_ptr[2][1]=4;
g->cell_line_node_ptr[2][2]=2;
g->cell_line_node_ptr[2][3]=1;
g->cell_line_node_ptr[2][4]=3;
g->cell_line_node_ptr[3][1]=2;
g->cell_line_node_ptr[3][2]=3;
g->cell_line_node_ptr[3][3]=4;
g->cell_line_node_ptr[3][4]=1;
g->cell_line_node_ptr[4][1]=1;
g->cell_line_node_ptr[4][2]=2;
g->cell_line_node_ptr[4][3]=3;
g->cell_line_node_ptr[4][4]=4;
g->cell_line_node_ptr[5][1]=1;
g->cell_line_node_ptr[5][2]=3;
g->cell_line_node_ptr[5][3]=2;
g->cell_line_node_ptr[5][4]=4;
g->cell_line_node_ptr[6][1]=1;
g->cell_line_node_ptr[6][2]=4;
g->cell_line_node_ptr[6][3]=2;
g->cell_line_node_ptr[6][4]=3;
g->cell_face_node_ptr[1][1]=2;
g->cell_face_node_ptr[1][2]=3;
g->cell_face_node_ptr[1][3]=4;
g->cell_face_node_ptr[1][4]=1;
g->cell_face_node_ptr[2][1]=1;
g->cell_face_node_ptr[2][2]=3;
g->cell_face_node_ptr[2][3]=4;
g->cell_face_node_ptr[2][4]=2;
g->cell_face_node_ptr[3][1]=1;
g->cell_face_node_ptr[3][2]=2;
g->cell_face_node_ptr[3][3]=4;
g->cell_face_node_ptr[3][4]=3;
g->cell_face_node_ptr[4][1]=1;
g->cell_face_node_ptr[4][2]=2;
g->cell_face_node_ptr[4][3]=3;
g->cell_face_node_ptr[4][4]=4;
break;
}
g->cell_node_num= g->spacedim+1;
///////////////////////////////////////
///// Calculate min/max and node numbers ?!
/// MESH*::noOfPoints is "private (at least in principle)" !!
g->prev_node_num=0;
g->xmin=g->ymin=g->zmin= 1.0e30;
g->xmax=g->ymax=g->zmax= -1.0e30;
switch (spacedim)
{
case 1:
{
PT1 *patch;
DListIter<PT1> iter(m->castToMESH1()->points());
while ((patch = iter.all()))
{
if (patch->x< g->xmin) g->xmin=patch->x;
if (patch->x> g->xmax) g->xmax=patch->x;
}
g->ymin=g->zmin=0.0;
g->ymax=g->zmax=0.0;
}
break;
case 2:
{
PT2 *patch;
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?