Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F77536496
mesher.c
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, Aug 15, 01:25
Size
11 KB
Mime Type
text/x-c
Expires
Sat, Aug 17, 01:25 (2 d)
Engine
blob
Format
Raw Data
Handle
19624441
Attached To
rAKA akantu
mesher.c
View Options
/*
Copyright 2008 Guillaume ANCIAUX (guillaume.anciaux@epfl.ch)
This file is part of ParaViewHelper.
ParaViewHelper is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
ParaViewHelper is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ParaViewHelper. If not, see <http://www.gnu.org/licenses/>.
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
static double * coordinates;
static int * connectivity;
inline int compareCoords(double * c1,double * c2){
double tol=0.000001;
double x1,x2,y1,y2,z1,z2;
x1 = c1[0];x2 = c2[0];
if (fabs(x1-x2) > tol && x1 < x2) return -1;
if (fabs(x1-x2) > tol && x2 < x1) return 1;
y1 = c1[1];y2 = c2[1];
if (fabs (y1-y2) > tol && y1 < y2) return -1;
if (fabs (y1-y2) > tol && y2 < y1) return 1;
z1 = c1[2];z2 = c2[2];
if (fabs (z1-z2) > tol && z1 < z2) return -1;
if (fabs (z1-z2) > tol && z2 < z1) return 1;
return 0;
}
inline void AddNode(double * coords,int * nb_node, double x,double y,double z){
coords[3*(*nb_node)+0]=x;
coords[3*(*nb_node)+1]=y;
coords[3*(*nb_node)+2]=z;
(*nb_node)++;
}
inline void AddMedianNode(double *coords,int * nb_node,int n1,int n2){
coords[3*(*nb_node)+0]= (coords[3*n1+0]+coords[3*n2+0])/2.;
coords[3*(*nb_node)+1]= (coords[3*n1+1]+coords[3*n2+1])/2.;
coords[3*(*nb_node)+2]= (coords[3*n1+2]+coords[3*n2+2])/2.;
(*nb_node)++;
}
inline void AddElement(int * countelements,int * indexes,
int n0,int n1,int n2,
int n3,int n4,int n5,int n6,int n7,int n8,int n9){
connectivity[10*(*countelements)+0]=indexes[n0]+1;
connectivity[10*(*countelements)+1]=indexes[n1]+1;
connectivity[10*(*countelements)+2]=indexes[n2]+1;
connectivity[10*(*countelements)+3]=indexes[n3]+1;
connectivity[10*(*countelements)+4]=indexes[n4]+1;
connectivity[10*(*countelements)+5]=indexes[n5]+1;
connectivity[10*(*countelements)+6]=indexes[n6]+1;
connectivity[10*(*countelements)+7]=indexes[n7]+1;
connectivity[10*(*countelements)+8]=indexes[n8]+1;
connectivity[10*(*countelements)+9]=indexes[n9]+1;
(*countelements)++;
}
inline int AddInSorted(double * node_coords,double * coords_array,int * index_array,int N){
int low = 0;
int high = N - 1;
int mid=0;
int i = 0;
if (N == 0){
coords_array[0] = node_coords[0];
coords_array[1] = node_coords[1];
coords_array[2] = node_coords[2];
index_array[0] = -1;
return 0;
}
while (low <= high)
{
mid = (low + high) / 2;
switch (compareCoords(node_coords,coords_array+3*mid)){
case -1: high = mid -1; break;
case 1: low = mid + 1; break;
case 0: return mid;break;//such a node as already been inserted
default: fprintf(stderr,"AAAAAAAAARRRRRRRRRRG !! should not append\n");
}
}
mid = (low + high) / 2;
if (compareCoords(coords_array+3*mid,node_coords)==-1)
++mid;
//shift all elements to insert at the right place the element we want
double * c1 = coords_array+3*mid;
double * c2;
int * i1 = index_array+mid;
int * i2;
for (i = N-1; i >= mid ; --i){
c1 = coords_array+3*i;
c2 = c1+3;
i1 = index_array+i;
i2 = i1+1;
c2[0] = c1[0];
c2[1] = c1[1];
c2[2] = c1[2];
*i2 = * i1;
}
c1[0] = node_coords[0];
c1[1] = node_coords[1];
c1[2] = node_coords[2];
*i1 = -1;
return mid;
}
inline void AddNodesInGlobal(double * coords,int size,int * global_size,double * sorted_coordinates,int * sorted_indexes,int * indexes){
// int size2 = *global_size;
int i;
double c[3];
for (i = 0 ; i < size ; ++i){
//speedup search in sorted list to see if node already exists
c[0] = coords[3*i];c[1] = coords[3*i+1];c[2] = coords[3*i+2];
// if (c[0] == .5 && c[1] == 1 && c[2] == 1) debug_flag = 1;
int ind = AddInSorted(c,sorted_coordinates,sorted_indexes,*global_size);
if (sorted_indexes[ind] == -1) {//node does not already exist
coordinates[3*(*global_size)] = coords[3*i];
coordinates[3*(*global_size)+1] = coords[3*i+1];
coordinates[3*(*global_size)+2] = coords[3*i+2];
indexes[i] = (*global_size);
sorted_indexes[ind] = (*global_size);
(*global_size)++;
}
else {//node already exist
indexes[i] = sorted_indexes[ind];
}
/* fprintf(stderr,"rnodes: start\n"); */
/* for (j = 0 ; j < *global_size ; ++j){ */
/* fprintf(stderr,"(%e,%e,%e,%d)[%d]\n", */
/* sorted_coordinates[3*j],sorted_coordinates[3*j+1],sorted_coordinates[3*j+2], */
/* sorted_indexes[j],j); */
/* } */
/* fprintf(stderr,"end\n"); */
}
}
void write_mesh_UNV(const char* outputfile,int nodes,int elements){
int i;
FILE *op = fopen(outputfile, "w");
//nodal values (not used at present time)
int exp_coord_sys_num = 0;
int disp_coord_sys_num = 0;
int color = 0;
//element values
int fe_descriptor_id = 111;
int phys_prop_tab_num = 2;
int mat_prop_tab_num = 1;
int n_nodes = 4;
//start of node block
fprintf(op," -1\n 2411\n");
for (i=0;i<nodes;i++){
fprintf(op," %d %d %d %d\n",
i,exp_coord_sys_num,disp_coord_sys_num,color);
fprintf(op," %.15e %.15e %.15e\n",
coordinates[3*i+0],
coordinates[3*i+1],
coordinates[3*i+2]);
}
//end of node block
fprintf(op," -1\n");
//start of element block
fprintf(op," -1\n 2412\n");
color = 7;
for (i=0;i<elements;i++){
fprintf(op," %d %d %d %d %d %d\n",
i,fe_descriptor_id,phys_prop_tab_num,mat_prop_tab_num,color,n_nodes);
fprintf(op," %d %d %d %d\n",
connectivity[10*i+0]-1,connectivity[10*i+1]-1,connectivity[10*i+2]-1,connectivity[10*i+3]-1);
}
//end of element block
fprintf(op,"-1\n");
fclose(op);
}
void BuildGridMesh(double ** pos,int ** conn,
int * nodes,int * elements,
double L, double D, double W,
double NL,double ND, double NW)
{
int i,j,k;
int countnodes,countelements;
int nstart,n0,n1,n2,n3,n4,n5,n6,n7,n8,n9;
double lcube,dcube,wcube;
double x0,y0,z0;
// double tol = 0.00001;
lcube = L/NL;
dcube = D/ND;
wcube = W/NW;
*nodes = 35*NL*ND*NW; /* 35 is the number of nodes per unit cube; node is overestimated but identical nodes will be removed later on */
*elements = 12*NL*ND*NW;
*pos = malloc(*nodes*3*sizeof(double));
coordinates = *pos;
double * sorted_coordinates = malloc(*nodes*3*sizeof(double));
int * sorted_indexes = malloc(*nodes*sizeof(int));
*conn = malloc(*elements*10*sizeof(int));
connectivity = *conn;
// int global_size=0;
countnodes=0;
countelements=0;
for (i=0;i<NW;i++)
for (j=0;j<ND;j++)
for (k=0;k<NL;k++)
{
if (countelements%1000 == 0)
fprintf(stderr,"building element %d/%d\n",countelements,*elements);
x0 = lcube*k;
y0 = dcube*j;
z0 = wcube*i;
double temp_coords[3*35];
int indexes[35];
int cpt=0;
/* Base corner nodes */
nstart = countnodes+1;
AddNode(temp_coords,&cpt,x0,y0,z0);
AddNode(temp_coords,&cpt,x0+lcube,y0,z0);
AddNode(temp_coords,&cpt,x0+lcube,y0+dcube,z0);
AddNode(temp_coords,&cpt,x0,y0+dcube,z0);
/* Top corner nodes */
AddNode(temp_coords,&cpt,x0,y0,z0+wcube);
AddNode(temp_coords,&cpt,x0+lcube,y0,z0+wcube);
AddNode(temp_coords,&cpt,x0+lcube,y0+dcube,z0+wcube);
AddNode(temp_coords,&cpt,x0,y0+dcube,z0+wcube);
/* Mid node */
n1=0;n2=6;
AddMedianNode(temp_coords,&cpt,n1,n2);
/* Mid nodes on the faces */
n1=0;n2=1;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=1;n2=2;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=2;n2=3;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=0;n2=3;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=4;n2=5;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=5;n2=6;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=6;n2=7;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=7;n2=4;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=0;n2=4;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=1;n2=5;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=2;n2=6;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=3;n2=7;
AddMedianNode(temp_coords,&cpt,n1,n2);
/* Mid nodes on the diagonals of the faces */
n1=4;n2=1;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=6;n2=1;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=7;n2=2;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=7;n2=0;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=3;n2=1;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=7;n2=5;
AddMedianNode(temp_coords,&cpt,n1,n2);
/* Joining the corner nodes to the central mid node n9 */
n1=0;n2=8;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=1;n2=8;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=2;n2=8;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=3;n2=8;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=4;n2=8;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=5;n2=8;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=6;n2=8;
AddMedianNode(temp_coords,&cpt,n1,n2);
n1=7;n2=8;
AddMedianNode(temp_coords,&cpt,n1,n2);
//memcpy(coordinates+(nstart-1)*3,temp_coords,3*nb_nodes*sizeof(double));
AddNodesInGlobal(temp_coords,cpt,&countnodes,sorted_coordinates,sorted_indexes,indexes);
// countnodes += cpt;
/* CONNECTIVITY */
n0=0;n1=4;n2=1;n3=8;n4=17;n5=21;n6=9;n7=27;n8=31;n9=28;
AddElement(&countelements,indexes,n0,n1,n2,n3,n4,n5,n6,n7,n8,n9);
n0=1;n1=4;n2=5;n3=8;n4=21;n5=13;n6=18;n7=28;n8=31;n9=32;
AddElement(&countelements,indexes,n0,n1,n2,n3,n4,n5,n6,n7,n8,n9);
n0=1;n1=5;n2=6;n3=8;n4=18;n5=14;n6=22;n7=28;n8=32;n9=33;
AddElement(&countelements,indexes,n0,n1,n2,n3,n4,n5,n6,n7,n8,n9);
n0=1;n1=6;n2=2;n3=8;n4=22;n5=19;n6=10;n7=28;n8=33;n9=29;
AddElement(&countelements,indexes,n0,n1,n2,n3,n4,n5,n6,n7,n8,n9);
n0=2;n1=6;n2=7;n3=8;n4=19;n5=15;n6=23;n7=29;n8=33;n9=34;
AddElement(&countelements,indexes,n0,n1,n2,n3,n4,n5,n6,n7,n8,n9);
n0=2;n1=7;n2=3;n3=8;n4=23;n5=20;n6=11;n7=29;n8=34;n9=30;
AddElement(&countelements,indexes,n0,n1,n2,n3,n4,n5,n6,n7,n8,n9);
n0=3;n1=7;n2=0;n3=8;n4=20;n5=24;n6=12;n7=30;n8=34;n9=27;
AddElement(&countelements,indexes,n0,n1,n2,n3,n4,n5,n6,n7,n8,n9);
n0=7;n1=4;n2=0;n3=8;n4=16;n5=17;n6=24;n7=34;n8=31;n9=27;
AddElement(&countelements,indexes,n0,n1,n2,n3,n4,n5,n6,n7,n8,n9);
n0=4;n1=7;n2=5;n3=8;n4=16;n5=26;n6=13;n7=31;n8=34;n9=32;
AddElement(&countelements,indexes,n0,n1,n2,n3,n4,n5,n6,n7,n8,n9);
n0=7;n1=6;n2=5;n3=8;n4=15;n5=14;n6=26;n7=34;n8=33;n9=32;
AddElement(&countelements,indexes,n0,n1,n2,n3,n4,n5,n6,n7,n8,n9);
n0=0;n1=1;n2=3;n3=8;n4=9;n5=25;n6=12;n7=27;n8=28;n9=30;
AddElement(&countelements,indexes,n0,n1,n2,n3,n4,n5,n6,n7,n8,n9);
n0=1;n1=2;n2=3;n3=8;n4=10;n5=11;n6=25;n7=28;n8=29;n9=30;
AddElement(&countelements,indexes,n0,n1,n2,n3,n4,n5,n6,n7,n8,n9);
}
*nodes = countnodes;
// write_mesh_UNV("cube.unv",*nodes,*elements);
}
Event Timeline
Log In to Comment