wzsxgrid.cxx
来自「有限元学习研究用源代码(老外的),供科研人员参考」· CXX 代码 · 共 1,773 行 · 第 1/3 页
CXX
1,773 行
DListIter<PT2> iter(m->castToMESH2()->points());
while ((patch = iter.all()))
{
if (patch->x< g->xmin) g->xmin=patch->x;
if (patch->x> g->xmax) g->xmax=patch->x;
if (patch->y< g->ymin) g->ymin=patch->y;
if (patch->y> g->ymax) g->ymax=patch->y;
}
g->zmin=0.0;
g->zmax=0.0;
}
break;
case 3:
{
PT3 *patch;
DListIter<PT3> iter(m->castToMESH3()->points());
while ((patch =iter.all()))
{
if (patch->x< g->xmin) g->xmin=patch->x;
if (patch->x> g->xmax) g->xmax=patch->x;
if (patch->y< g->ymin) g->ymin=patch->y;
if (patch->y> g->ymax) g->ymax=patch->y;
if (patch->z< g->zmin) g->zmin=patch->z;
if (patch->z> g->zmax) g->zmax=patch->z;
}
}
break;
}
///////////////////////////////////////
///// parse cell materials
PATCH *patch;
i=0;
g->cell_mat_num=0;
m->resetElemIter();
while ((patch = m->elemIterAll()))
{
if (patch->Class()>g->cell_mat_num) g->cell_mat_num=patch->Class();
}
g->cell_mat_cell_num=(int*)calloc(g->cell_mat_num+1,sizeof(int));
/////////////////////////////////////////////////////
///// parse face materials (boundary condition types)
g->face_node_num=g->spacedim;
g->face_mat_num=0;
switch(spacedim)
{
case 1:
{
PT1* patch;
DListIter<PT1> iter(m->castToMESH1()->points());
while ((patch = iter.all()))
if (patch->onBoundary()&&patch->Class())
{
if (patch->Class()>g->face_mat_num)
g->face_mat_num=patch->Class();
}
}
break;
case 2:
{
EDG2* patch;
DListIter<EDG2> iter(m->castToMESH2()->lines());
i=0;
while ((patch = iter.all()))
if (patch->onBoundary()&&patch->Class())
{
if (patch->Class()>g->face_mat_num)
g->face_mat_num=patch->Class();
}
}
break;
case 3:
{
TR3* patch;
DListIter<TR3> iter(m->castToMESH3()->triangles());;
while ((patch = iter.all()))
if (patch->onBoundary()&&patch->Class())
{
if (patch->Class()>g->face_mat_num)
g->face_mat_num=patch->Class();
}
break;
}
}
g->face_mat_face_num=(int*) calloc(g->face_mat_num+1,sizeof(int));
number_stuff(g);
sxgPrintData(g);
dprintf(inf)("sxGrid created from mesh\n",meshname);
return g;
}
void sxDestroy(sxGrid g)
{
dprintf(inf)("Deleting sxgrid...\n");
free(g->cell_mat_cell_num);
free(g->face_mat_face_num);
delete g->mesh;
delete g;
dprintf(inf)("Deletion done.\n");
}
/////////////////////////////////////////
// Interface to mesh refinement
/////////////////////////////////////////
// methods on simplices
int sxDepth(Simplex sx)
{
// My depth is relative to the finest level!!!
return sx->grid->mesh->maxDepth-sx->priv->patch->Depth();
}
void sxMark(Simplex sx)
{
sx->priv->patch->setMark();
}
void sxUnMark(Simplex sx)
{
sx->priv->patch->unMark();
}
int sxMarked(Simplex sx)
{
if (sx->priv->patch->marked()) return 1; else return 0;
}
//////// methods on mesh
int sxMaxDepth(sxGrid g)
{
return g->mesh->maxDepth;
}
void sxRefine(sxGrid g)
{
MESH *m=g->mesh;
int i;
dprintf(inf)("Refining sxGrid...\n");
m->Refine();
g->prev_node_num=g->node_num;
number_stuff(g);
if (debug(inf)) sxgPrintData(g);
dprintf(inf)("Refinement done\n");
}
void sxFlatten(sxGrid g)
{
g->mesh->Flatten();
}
//// Interpolation of vectors
/////////////////////////////////////////
static double gshape[30];
static int sx_nodes[5]={1,2,3,4};
#define COORD(DIM,NODE,DIRm1) coord[1+DIRm1+(cellnodes[NODE]-1)*DIM]
void sxFVMFactors(Simplex sx)
{
sx->priv->FVMFactors(
sx->priv->coord,
sx_nodes,
&sx->volume,
sx->nodefac+1,
sx->linefac+1
);
}
void sxFVMMatrices(Simplex sx)
{
register int i,j;
sx->priv->FVMFactors(
sx->priv->coord,
sx_nodes,
&sx->volume,
sx->nodefac+1,
sx->linefac+1
);
for(i=1;i<=sx->nnodes;i++)
for(j=1;j<=sx->nnodes;j++)
sx->stiffmat[i][j]=sx->massmat[i][j]=0.0;
for(i=1;i<=sx->nlines;i++)
{
sx->stiffmat[sx->linenodes[i][1]][sx->linenodes[i][2]]=-sx->linefac[i];
sx->stiffmat[sx->linenodes[i][2]][sx->linenodes[i][1]]=-sx->linefac[i];
sx->stiffmat[sx->linenodes[i][1]][sx->linenodes[i][1]]+=sx->linefac[i];
sx->stiffmat[sx->linenodes[i][2]][sx->linenodes[i][2]]+=sx->linefac[i];
}
for(i=1;i<=sx->nnodes;i++)
sx->massmat[i][i]=sx->nodefac[i];
}
void sxFEMFactors(Simplex sx)
{
register int id,in,ie;
sx->priv->FEMFactors(
sx->priv->coord,
sx_nodes,
&sx->volume,
sx->nodefac+1,
sx->linefac+1,
gshape
);
for(ie=1;ie<=sx->nlines;ie++)
sx->linefac[ie]= -sx->linefac[ie];
for(in=1;in<=sx->nnodes;in++)
sx->nodefac[in]=sx->volume/((double)sx->nnodes);
for(in=1;in<=sx->nnodes;in++)
for(id=1;id<=sx->spacedim;id++)
sx->dgrad[in][id]=gshape[(in-1)*sx->spacedim+(id-1)];
}
static double xn[4]=
{
0,
0.1666666666,
0.0833333333,
0.05
};
void sxFEMMatrices(Simplex sx)
{
register int i,j,in,id;
double xv;
sx->priv->FVMFactors(
sx->priv->coord,
sx_nodes,
&sx->volume,
sx->nodefac+1,
sx->linefac+1
);
for(i=1;i<=sx->nnodes;i++)
sx->stiffmat[i][i]=0.0;
for(i=1;i<=sx->nlines;i++)
{
sx->stiffmat[sx->linenodes[i][1]][sx->linenodes[i][2]]=-sx->linefac[i];
sx->stiffmat[sx->linenodes[i][2]][sx->linenodes[i][1]]=-sx->linefac[i];
sx->stiffmat[sx->linenodes[i][1]][sx->linenodes[i][1]]+=sx->linefac[i];
sx->stiffmat[sx->linenodes[i][2]][sx->linenodes[i][2]]+=sx->linefac[i];
}
/* npar is wrong in celltye!!!*/
xv=xn[sx->dim]*sx->volume;
for(i=1;i<=sx->nnodes;i++)
for(j=1;j<=sx->nnodes;j++)
sx->massmat[i][j]=xv;
for(i=1;i<=sx->nnodes;i++)
sx->massmat[i][i]+=xv;
for(in=1;in<=sx->nnodes;in++)
for(id=1;id<=sx->dim;id++)
sx->dgrad[in][id]=gshape[(in-1)*sx->spacedim+(id-1)];
}
void sxFVEMFactors(Simplex sx)
{
register int id,in;
sx->priv->FEMFactors(
sx->priv->coord,
sx_nodes,
&sx->volume,
sx->nodefac+1,
sx->linefac+1,
gshape
);
for(in=1;in<=sx->nnodes;in++)
for(id=1;id<=sx->spacedim;id++)
sx->dgrad[in][id]=gshape[(in-1)*sx->spacedim+(id-1)];
sx->priv->FVMFactors(
sx->priv->coord,
sx_nodes,
&sx->volume,
sx->nodefac+1,
sx->linefac+1
);
}
double* sxCenter(Simplex sx, double *ccrd)
{
register int i,idim;
for (idim=1;idim<=sx->spacedim;idim++) ccrd[idim]=0.0;
for(i=1;i<=sx->nnodes;i++)
{
for (idim=1;idim<=sx->spacedim;idim++)
ccrd[idim]+=sx->coord[i][idim];
}
for (idim=1;idim<=sx->spacedim;idim++)
ccrd[idim]/=(double)(sx->nnodes);
return ccrd;
}
////////////////////////////////////////////////////////////////
sxFunction sxCreateFunction(sxGrid g, char *name,int num, int comp)
{
double *d=NULL;
switch (comp)
{
case SX_NODES:
d=new double[num*g->node_num+1]; return d;
case SX_CELLS:
d=new double[num*g->cell_num+1]; return d;
case SX_EDGES:
d=new double[num*g->line_num+1]; return d;
return NULL;
case SX_ALL_FACES:
d=new double[num*g->all_face_num+1]; return d;
return NULL;
default:
dprintf(err)("Function type %d not supported",comp);
}
return NULL;
}
void sxDestroyFunction(sxGrid g, sxFunction f)
{
delete f;
}
// this interpolates _only_ from old level to new level
sxFunction sxInterpolateFunctionFromPrevLevel
(sxGrid g, sxFunction f_old, int num, int comp)
{
double *f_new=NULL;
MESH *mesh=g->mesh;
if (num>1)
{
dprintf(err)("sorry, interpolation for multiple functions currently not implemented\n");
return NULL;
}
switch (comp)
{
case SX_NODES:
int i;
// inject coarse mesh values
// assume old values come first!?
f_new=new double[num*g->node_num+1];
for (i=1;i<=g->prev_node_num;i++) f_new[i]=f_old[i];
for (;i<=g->node_num;i++) f_new[i]=0.0;
// MLNode seems to be undefined unless
// LinElemInt:: setMLNodeNumbers has been called!
// {
// int i, node_old,node_new;
// PATCH* p;
// f_new=new double[num*g->node_num+1];
// for (i=1;i<=g->node_num;i++) f_new[i]=0.0;
//
// mesh->resetPtIter();
// while (p=mesh->ptIterAll())
// if (p->Depth()<mesh->maxDepth)
// {
// node_new = p->getNode();
// node_old = p->MLNode(1);
// f_new[node_new] = f_old[node_old];
// }
// }
//interplolate new values
{
PATCH* ed;
PATCH *pm;
int node1, node2, pmNode;
int iint=0;
mesh->resetEdgIter();
while ((ed=mesh->edgIterAllOfHistory()))
{
if (ed->refined())
{
pm = ed->getMidPoint();
pmNode = pm->getNode();
if (pmNode>g->prev_node_num)
{
iint++;
const Vector<PT*>& p=ed->getPoints();
node1 = p[1]->getNode();
node2 = p[2]->getNode();
f_new[pmNode] = 0.5*(f_new[node1] + f_new[node2]);
}
}
}
}
delete f_old;
return f_new;
break;
case SX_CELLS:
dprintf(err)("sorry, interpolation of cell functions currently not implemented\n");
return NULL;
break;
case SX_EDGES:
dprintf(err)("sorry, interpolation of line functions currently not implemented\n");
return NULL;
break;
case SX_ALL_FACES:
case SX_FACES:
dprintf(err)("sorry, interpolation of face functions currently not implemented\n");
return NULL;
}
return NULL;
}
int sxAdaptive(sxGrid g)
{
return 1;
}
//get data of cell with x in it.
//void sxFindCell(SimplexGrid g,
// double coord[],
// double baryc_coord[],
// SimplexData *sx)
//{
//
// p=findPatch(coord,baryc_coord);
// SetCellData(sx,p);
// baryc_coord[2]=1.0-baryc_coord[1]-baryc_coord[0];
// return;
//
//}
////////////////////////////////////////////////////////////////
enum pauseMode Pause()
{
static Bool Continue = False;
if (Continue == True) return noPicture;
if (!Cmd.isTrue("pause")) { Continue = True; return noPicture; }
char s[5]; cout << " <CR>"; cout.flush(); gets(s); strToLower(s);
if (strchr(s,'q') || strchr(s,'e')) { cout << "\nEXIT FORCED\n"; exit(1); }
if (strchr(s,'c') || strchr(s,'g')) Continue = True;
return noPicture;
}
void sxLocatePoint(sxGrid g, double *coord, Simplex sx, double *bcoord)
{
Vector<Real>x(4); Vector<Real>xUnit(4);
for(int i=1;i<=g->spacedim;i++)
x[i]=coord[i];
switch (g->spacedim)
{
case 1:
EDG1 * cell1;
cell1 = (EDG1*) g->mesh->findPatch(x,xUnit,cell1);
if (cell1==0) {sx->index=0; return;}
setCell1Data(sx,cell1);
bcoord[1]=1.0-xUnit[1];
bcoord[2]=xUnit[1];
break;
case 2:
TR2 * cell2;
cell2=(TR2*)g->mesh->findPatch(x,xUnit,cell2);
if (cell2==0) {sx->index=0; return;}
// this order is somehow hard coded, but how
// at least we can look it up in elements2.cc!
bcoord[1]=1.0-xUnit[1]-xUnit[2];
bcoord[2]=xUnit[1];
bcoord[3]=xUnit[2];
setCell2Data(sx,cell2);
break;
case 3:
TET3 * cell3;
cell3=(TET3*)g->mesh->findPatch(x,xUnit,cell3);
if (cell3==0) {sx->index=0; return;}
setCell3Data(sx,cell3);
bcoord[1]=1.0-xUnit[1]-xUnit[2]-xUnit[3];
bcoord[2]=xUnit[1];
bcoord[3]=xUnit[2];
bcoord[4]=xUnit[3];
break;
}
}
sxFunction sxInterpolateFunctionFromPrevGrid
(sxGrid g, sxGrid g_old, sxFunction f_old, int num, int comp)
{
double *f_new=NULL;
MESH *mesh=g->mesh;
if (num>1)
{
dprintf(err)("sorry, interpolation for multiple functions currently not implemented\n");
return NULL;
}
switch (comp)
{
case SX_NODES:
int i;
SimplexData sx,sx_old;
sx_old.priv=new PrivateSimplexData;
for(i=1;i<=g->spacedim+1;i++)
sx_old.coord[i]=sx_old.priv->coord+(i-1)*g->spacedim;
sx_old.node= new int[5];
sx_old.line= new int[5];
double bcoord[5];
f_new=new double[num*g->node_num+1];
SX_LOOP(g,sx,SX_NODES)
{
sxLocatePoint(g_old,sx.coord[1],&sx_old,bcoord);
double f=0.0;
for (int in=1;in<=g->spacedim+1;in++)
// cout<< bcoord[in]<< ' ' << f_old[sx_old.node[in]]<<endl<<flush,
f+=bcoord[in]*f_old[sx_old.node[in]];
// cout<<endl;
f_new[sx.index]=f;
}
delete[] sx_old.node;
delete[] sx_old.line;
delete sx_old.priv;
return f_new;
break;
case SX_CELLS:
dprintf(err)("sorry, interpolation of cell functions currently not implemented\n");
return NULL;
break;
case SX_EDGES:
dprintf(err)("sorry, interpolation of line functions currently not implemented\n");
return NULL;
break;
case SX_ALL_FACES:
case SX_FACES:
dprintf(err)("sorry, interpolation of face functions currently not implemented\n");
return NULL;
}
return NULL;
}
sxGrid global_sxg=NULL;
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?