59 bool read_binary =
fileExists(basename +
".bpts");
60 std::string pts_file = read_binary ? basename +
".bpts" : basename +
".pts";
62 FILE* pts_fd = fopen(pts_file.c_str(),
"r");
64 fprintf(stderr,
"Error: could not open file: %s. Aborting!\n", pts_file.c_str());
70 if(read_binary) fread(buffer,
sizeof(
char),
HDR_SIZE, pts_fd);
71 else fgets( buffer, 2048, pts_fd);
72 sscanf(buffer,
"%lu", &numpts);
90 inline void read_headers(FILE* ele_fd, FILE* fib_fd,
bool read_binary,
size_t & numelem,
bool & twoFib)
97 fread(buffer,
sizeof(
char),
HDR_SIZE, ele_fd);
102 sscanf(buffer,
"%lu", &numelem);
107 fread(buffer,
sizeof(
char),
HDR_SIZE, fib_fd);
108 sscanf(buffer,
"%d", &numaxes);
114 short nread = sscanf(buffer,
"%f %f %f %f %f %f", fbuff, fbuff+1, fbuff+2, fbuff+3, fbuff+4, fbuff+5);
116 case 1: numaxes = fbuff[0];
break;
117 case 3: numaxes = 1;
break;
118 case 6: numaxes = 2;
break;
120 fprintf(stderr,
"%s error: read %d values in line 1 of fiber file." 121 "This is unexpected! Setting number of fibers to 1.\n",
129 fseek(fib_fd, 0, SEEK_SET);
132 twoFib = numaxes == 2;
147 inline void write_elem_headers(FILE* & ele_fd, FILE* & fib_fd,
bool binary,
size_t numelem,
bool twoFib)
153 int numfibres = twoFib ? 2 : 1;
158 snprintf(header,
sizeof header,
"%lu %d %d", numelem,
MT_ENDIANNESS, checksum);
159 fwrite(header,
HDR_SIZE,
sizeof(
char), ele_fd);
162 fprintf(ele_fd,
"%lu\n", numelem);
167 snprintf(header,
sizeof header,
"%d %lu %d %d", numfibres, (
unsigned long int)numelem,
MT_ENDIANNESS, checksum);
168 fwrite(header,
sizeof(
char),
HDR_SIZE, fib_fd);
171 fprintf(fib_fd,
"%d\n", numfibres);
192 snprintf(header,
sizeof header,
"%lu %d %d", numpts,
MT_ENDIANNESS, checksum);
193 fwrite(header,
HDR_SIZE,
sizeof(
char), pts_fd);
196 fprintf(pts_fd,
"%lu\n", numpts);
209 template<
class T,
class S>
212 const int bufsize = 2048;
213 const int max_elem_size = 9;
215 char buffer[bufsize];
217 int n[max_elem_size];
222 mesh.
dsp.resize(bsize+1);
223 mesh.
tag.resize(bsize);
225 mesh.
con.resize(bsize*max_elem_size);
227 size_t ele_read = 0, con_size = 0;
229 T* elem = mesh.
con.data();
231 for(
size_t i=0; i<bsize; i++)
235 size_t r = fread(&etbuff,
sizeof(
int), 1, fd);
240 char* ptr = fgets( buffer, bufsize, fd);
241 if(ptr == NULL)
break;
242 sscanf(buffer,
"%s %d %d %d %d %d %d %d %d %d", etype_str, n, n+1, n+2, n+3, n+4, n+5, n+6, n+7, n+8);
278 fprintf(stderr,
"Error: Unsupported element type!\n");
281 mesh.
dsp[i+1] = mesh.
dsp[i] + nodes;
286 fread(n,
sizeof(
int), nodes+1, fd);
288 for(
int j=0; j<nodes; j++) elem[j] = n[j];
290 mesh.
tag[i] = n[nodes];
291 ref_eidx[i] = bstart + i;
296 mesh.
dsp.resize(ele_read+1);
297 mesh.
tag.resize(ele_read);
299 ref_eidx.
resize(ele_read);
300 mesh.
con.resize(con_size);
313 template<
class T,
class S>
316 const T* con = mesh.
con.data();
322 for(
size_t eidx=0; eidx < mesh.
l_numelem; eidx++)
324 fwrite(&mesh.
type[eidx], 1,
sizeof(
elem_t), fd);
325 T esize = mesh.
dsp[eidx+1] - mesh.
dsp[eidx];
327 wbuff[esize] = mesh.
tag[eidx];
329 fwrite(wbuff, esize+1,
sizeof(
int), fd);
335 for(
size_t eidx=0; eidx < mesh.
l_numelem; eidx++)
337 switch(mesh.
type[eidx])
340 fprintf(fd,
"Ln %d %d %d\n", con[0], con[1], mesh.
tag[eidx]);
345 fprintf(fd,
"Tr %d %d %d %d\n", con[0], con[1], con[2], mesh.
tag[eidx]);
350 fprintf(fd,
"Qd %d %d %d %d %d\n", con[0], con[1], con[2], con[3], mesh.
tag[eidx]);
355 fprintf(fd,
"Tt %d %d %d %d %d\n", con[0], con[1], con[2], con[3], mesh.
tag[eidx]);
360 fprintf(fd,
"Py %d %d %d %d %d %d\n", con[0], con[1], con[2], con[3],
361 con[4], mesh.
tag[eidx]);
366 fprintf(fd,
"Pr %d %d %d %d %d %d %d\n", con[0], con[1], con[2], con[3],
367 con[4], con[5], mesh.
tag[eidx]);
372 fprintf(fd,
"Hx %d %d %d %d %d %d %d %d %d\n", con[0], con[1], con[2], con[3],
373 con[4], con[5], con[6], con[7], mesh.
tag[eidx]);
378 fprintf(stderr,
"Error: Unsupported element type!\n");
395 template<
class T,
class S>
398 const int bufsize = 2048;
399 char buffer[bufsize];
408 for(
size_t i=0; i<bsize; i++) {
410 size_t r = fread(fib,
sizeof(
float), 6, fd);
414 char* ptr = fgets( buffer, bufsize, fd);
415 if(ptr == NULL)
break;
416 sscanf(buffer,
"%f %f %f %f %f %f", fib, fib+1, fib+2, fib+3, fib+4, fib+5);
418 mesh.
fib[nr*3+0] = fib[0];
419 mesh.
fib[nr*3+1] = fib[1];
420 mesh.
fib[nr*3+2] = fib[2];
421 mesh.
she[nr*3+0] = fib[3];
422 mesh.
she[nr*3+1] = fib[4];
423 mesh.
she[nr*3+2] = fib[5];
436 for(
size_t i=0; i<bsize; i++) {
438 size_t r = fread(fib,
sizeof(
float), 3, fd);
442 char* ptr = fgets( buffer, bufsize, fd);
443 if(ptr == NULL)
break;
444 sscanf(buffer,
"%f %f %f", fib, fib+1, fib+2);
446 mesh.
fib[nr*3+0] = fib[0];
447 mesh.
fib[nr*3+1] = fib[1];
448 mesh.
fib[nr*3+2] = fib[2];
464 template<
class T,
class S>
474 fib[0] = mesh.
fib[i*3+0];
475 fib[1] = mesh.
fib[i*3+1];
476 fib[2] = mesh.
fib[i*3+2];
477 fib[3] = mesh.
she[i*3+0];
478 fib[4] = mesh.
she[i*3+1];
479 fib[5] = mesh.
she[i*3+2];
482 fwrite(fib, 6,
sizeof(
float), fd);
484 fprintf(fd,
"%f %f %f %f %f %f\n", fib[0], fib[1], fib[2], fib[3], fib[4], fib[5]);
491 fib[0] = mesh.
fib[i*3+0];
492 fib[1] = mesh.
fib[i*3+1];
493 fib[2] = mesh.
fib[i*3+2];
496 fwrite(fib, 3,
sizeof(
float), fd);
498 fprintf(fd,
"%f %f %f\n", fib[0], fib[1], fib[2]);
512 template<
class T,
class S>
516 MPI_Comm comm = mesh.
comm;
519 MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
521 FILE* ele_fd = NULL, *fib_fd = NULL;
522 size_t gnumelems = 0;
524 bool read_binary =
false;
529 read_binary =
fileExists(basename +
".belem");
530 std::string ele_file = read_binary ? basename +
".belem" : basename +
".elem";
531 std::string fib_file = read_binary ? basename +
".blon" : basename +
".lon";
533 ele_fd = fopen(ele_file.c_str(),
"r");
535 fprintf(stderr,
"Error: could not open file: %s. Aborting!\n", ele_file.c_str());
539 fib_fd = fopen(fib_file.c_str(),
"r");
541 fprintf(stderr,
"Error: could not open file: %s. Aborting!\n", fib_file.c_str());
546 MPI_Allreduce(MPI_IN_PLACE, &err, 1, MPI_INT, MPI_SUM, comm);
552 read_headers(ele_fd, fib_fd, read_binary, gnumelems, twoFib);
554 MPI_Bcast(&gnumelems,
sizeof(
size_t), MPI_BYTE, 0, comm);
555 MPI_Bcast(&twoFib,
sizeof(
bool), MPI_BYTE, 0, comm);
559 size_t blocksize = (gnumelems + size - 1) / size;
569 for(
int pid=1; pid<size; pid++) {
571 read_elem_block(ele_fd, read_binary, pid*blocksize, blocksize, meshbuff);
572 read_fib_block (fib_fd, read_binary, twoFib, blocksize, meshbuff);
577 MPI_Send(meshbuff.
dsp.data(), meshbuff.
dsp.size()*
sizeof(T), MPI_BYTE, pid,
SF_MPITAG, comm);
578 MPI_Send(meshbuff.
tag.data(), meshbuff.
tag.size()*
sizeof(T), MPI_BYTE, pid,
SF_MPITAG, comm);
579 MPI_Send(ref_eidx.
data(), ref_eidx.
size()*
sizeof(T), MPI_BYTE, pid,
SF_MPITAG, comm);
585 size_t con_size = meshbuff.
con.size();
586 MPI_Send(&con_size,
sizeof(
size_t), MPI_BYTE, pid,
SF_MPITAG, comm);
587 MPI_Send(meshbuff.
con.data(), con_size*
sizeof(T), MPI_BYTE, pid,
SF_MPITAG, comm);
608 MPI_Recv(mesh.
dsp.data(), mesh.
dsp.size()*
sizeof(T), MPI_BYTE, 0,
SF_MPITAG, comm, &stat);
609 MPI_Recv(mesh.
tag.data(), mesh.
tag.size()*
sizeof(T), MPI_BYTE, 0,
SF_MPITAG, comm, &stat);
610 MPI_Recv(ref_eidx.
data(), ref_eidx.
size()*
sizeof(T), MPI_BYTE, 0,
SF_MPITAG, comm, &stat);
617 MPI_Recv(&con_size,
sizeof(
size_t), MPI_BYTE, 0,
SF_MPITAG, comm, &stat);
618 mesh.
con.resize(con_size);
619 MPI_Recv(mesh.
con.data(), con_size*
sizeof(T), MPI_BYTE, 0,
SF_MPITAG, comm, &stat);
636 template<
class T,
class S>
639 const MPI_Comm comm = mesh.
comm;
642 MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
644 FILE* ele_fd = NULL, *fib_fd = NULL;
646 std::string ele_file = binary ? basename +
".belem" : basename +
".elem";
647 std::string fib_file = binary ? basename +
".blon" : basename +
".lon";
651 ele_fd = fopen(ele_file.c_str(),
"w");
653 fprintf(stderr,
"Error: could not open file: %s. Aborting!\n", ele_file.c_str());
656 fib_fd = fopen(fib_file.c_str(),
"w");
658 fprintf(stderr,
"Error: could not open file: %s. Aborting!\n", fib_file.c_str());
668 for(
int pid=0; pid < size; pid++)
672 ele_fd = fopen(ele_file.c_str(),
"a");
674 fprintf(stderr,
"Error: could not open file: %s. Aborting!\n", ele_file.c_str());
677 fib_fd = fopen(fib_file.c_str(),
"a");
679 fprintf(stderr,
"Error: could not open file: %s. Aborting!\n", fib_file.c_str());
694 template<
class T,
class S>
697 const MPI_Comm comm = surfmesh.
comm;
698 const bool binary =
false;
701 MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
707 ele_fd = fopen(surffile.c_str(),
"w");
709 fprintf(stderr,
"Error: could not open file: %s. Aborting!\n", surffile.c_str());
712 fprintf(ele_fd,
"%lu\n", surfmesh.
g_numelem);
717 for(
int pid=0; pid < size; pid++)
720 ele_fd = fopen(surffile.c_str(),
"a");
722 fprintf(stderr,
"Error: could not open file: %s. Aborting!\n", surffile.c_str());
745 const int bufsize = 2048;
746 char buffer[bufsize];
752 for(
size_t i=0; i<bsize; i++) {
754 size_t r = fread(pts,
sizeof(
float), 3, fd);
758 char* ptr = fgets( buffer, bufsize, fd);
759 if(ptr == NULL)
break;
760 sscanf(buffer,
"%f %f %f", pts, pts+1, pts+2);
762 xyz[nr*3+0] = pts[0];
763 xyz[nr*3+1] = pts[1];
764 xyz[nr*3+2] = pts[2];
782 size_t nnodes = xyz.
size() / 3;
785 for(
size_t i=0; i<nnodes; i++)
792 fwrite(pt, 3,
sizeof(
float), fd);
794 fprintf(fd,
"%f %f %f\n", pt[0], pt[1], pt[2]);
806 template<
class T,
class S>
810 MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
814 bool read_binary =
false;
819 std::string pts_file = read_binary ? basename +
".bpts" : basename +
".pts";
821 pts_fd = fopen(pts_file.c_str(),
"r");
823 fprintf(stderr,
"Error: could not open file: %s. Aborting!\n", pts_file.c_str());
829 if(read_binary) fread(buffer,
sizeof(
char),
HDR_SIZE, pts_fd);
830 else fgets( buffer, 2048, pts_fd);
831 sscanf(buffer,
"%lu", &gnumpts);
833 MPI_Bcast(&gnumpts,
sizeof(
size_t), MPI_BYTE, 0, comm);
836 size_t blocksize = (gnumpts + size - 1) / size;
844 for(
int pid = 1; pid < size; pid++) {
846 long int numsend = buff.
size();
848 MPI_Send(&numsend, 1, MPI_LONG, pid,
SF_MPITAG, comm);
849 MPI_Send(buff.
data(), numsend*
sizeof(S), MPI_BYTE, pid,
SF_MPITAG, comm);
854 long int numrecv = 0;
855 MPI_Recv(&numrecv, 1, MPI_LONG, 0,
SF_MPITAG, comm, &stat);
858 MPI_Recv(pts.
data(), numrecv*
sizeof(S), MPI_BYTE, 0,
SF_MPITAG, comm, &stat);
861 long int mysize = pts.
size() / 3;
865 interval(ptsidx, layout[rank], layout[rank+1]);
867 if(rank == 0) fclose(pts_fd);
877 template<
class T,
class S>
882 assert(meshlist.size() > 0);
883 assert(pts.
size() == (ptsidx.
size() * 3));
887 MPI_Comm comm = (*meshlist.begin())->comm;
890 MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
895 for(
int pid = 0; pid < size; pid++) {
896 size_t numsend = pts.
size();
897 MPI_Bcast(&numsend,
sizeof(
size_t), MPI_BYTE, pid, comm);
905 idxbuff.
resize(numsend / 3);
908 MPI_Bcast(ptsbuff.
data(), ptsbuff.
size()*
sizeof(S), MPI_BYTE, pid, comm);
909 MPI_Bcast(idxbuff.
data(), idxbuff.
size()*
sizeof(T), MPI_BYTE, pid, comm);
913 for(
auto it = meshlist.begin(); it != meshlist.end(); ++it)
921 for(
size_t i=0; i<rnod.
size(); i++)
929 for(
size_t j=0; j<idxbuff.
size(); j++) {
930 if(g2l.
count(idxbuff[j])) {
937 for(
size_t j=0; j<ridx.
size(); j++) {
938 T w = widx[j], r = ridx[j];
939 mesh.
xyz[w*3+0] = ptsbuff[r*3+0];
940 mesh.
xyz[w*3+1] = ptsbuff[r*3+1];
941 mesh.
xyz[w*3+2] = ptsbuff[r*3+2];
953 template<
class T,
class S>
956 FILE* vtk_file = fopen(file.c_str(),
"w");
957 if(vtk_file == NULL)
return;
960 fprintf (vtk_file,
"# vtk DataFile Version 3.0\n");
961 fprintf (vtk_file,
"vtk output\n");
962 fprintf (vtk_file,
"binary\n");
963 fprintf (vtk_file,
"DATASET UNSTRUCTURED_GRID\n\n");
964 fprintf (vtk_file,
"POINTS %lu float\n", mesh.
l_numpts);
969 for (
unsigned long int i=0; i<mesh.
l_numpts; i++ ) {
970 pts[0] =
htobe(p[0]);
971 pts[1] =
htobe(p[1]);
972 pts[2] =
htobe(p[2]);
973 fwrite(pts,
sizeof(
float), 3, vtk_file);
976 fprintf(vtk_file,
"\n");
978 fprintf (vtk_file,
"CELL_TYPES %lu\n", mesh.
l_numelem);
979 unsigned long int valcount = 0;
981 for(
unsigned long int i=0; i< mesh.
l_numelem; i++) {
1008 vtk_type =
htobe(vtk_type);
1009 fwrite(&vtk_type,
sizeof(
int), 1, vtk_file);
1010 valcount += mesh.
dsp[i+1] - mesh.
dsp[i] + 1;
1012 fprintf(vtk_file,
"\n");
1014 fprintf(vtk_file,
"CELLS %lu %lu\n", mesh.
l_numelem, valcount);
1015 const T* elem = mesh.
con.data();
1016 for(
unsigned long int i=0; i<mesh.
l_numelem; i++)
1018 int nodes = mesh.
dsp[i+1] - mesh.
dsp[i], n;
1019 int be_nodes =
htobe(nodes);
1020 fwrite(&be_nodes,
sizeof(
int), 1, vtk_file);
1022 for(
int j=0; j<nodes; j++) {
1023 n = elem[j]; n =
htobe(n);
1024 fwrite(&n,
sizeof(
int), 1, vtk_file);
1028 fprintf(vtk_file,
"\n");
1030 fprintf (vtk_file,
"CELL_DATA %lu \n", mesh.
l_numelem);
1031 fprintf (vtk_file,
"SCALARS elemTag int 1\n");
1032 fprintf (vtk_file,
"LOOKUP_TABLE default\n");
1033 for (
unsigned long int i=0; i<mesh.
l_numelem; i++ ) {
1035 fwrite(&t,
sizeof(
int), 1, vtk_file);
1041 fprintf (vtk_file,
"VECTORS fiber float\n");
1042 for (
unsigned long int i=0; i<mesh.
l_numelem; i++ ) {
1043 pts[0] =
htobe(
float(mesh.
fib[i*3+0])),
1044 pts[1] =
htobe(
float(mesh.
fib[i*3+1])),
1045 pts[2] =
htobe(
float(mesh.
fib[i*3+2]));
1046 fwrite(pts,
sizeof(
float), 3, vtk_file);
1048 fprintf(vtk_file,
"\n");
1051 fprintf (vtk_file,
"VECTORS sheet float\n");
1052 for (
unsigned long int i=0; i<mesh.
l_numelem; i++ ) {
1053 pts[0] =
htobe(
float(mesh.
she[i*3+0])),
1054 pts[1] =
htobe(
float(mesh.
she[i*3+1])),
1055 pts[2] =
htobe(
float(mesh.
she[i*3+2]));
1056 fwrite(pts,
sizeof(
float), 3, vtk_file);
1058 fprintf(vtk_file,
"\n");
The mesh storage class. It contains both element and vertex data.
The vector class and related algorithms.
void write_elements(const meshdata< T, S > &mesh, bool binary, std::string basename)
Read the element data (elements and fibers) of a CARP mesh.
void writeVTKmesh_binary(const meshdata< T, S > &mesh, std::string file)
Write a mesh in binary vtk format.
vector< T > & get_numbering(SF_nbr nbr_type)
Get the vector defining a certain numbering.
void write_elem_headers(FILE *&ele_fd, FILE *&fib_fd, bool binary, size_t numelem, bool twoFib)
Write the header of the element and fiber files.
void write_elem_block(FILE *fd, bool write_binary, const meshdata< T, S > &mesh)
Write the local element block to a file.
MPI_Comm comm
the parallel mesh is defined on a MPI world
void layout_from_count(const T count, vector< T > &layout, MPI_Comm comm)
void write_pts_block(FILE *&fd, bool write_binary, const vector< S > &xyz)
Write a chunk of points to a file.
vector< T > tag
element tag
vector< T > dsp
connectivity starting index of each element
void interval(vector< T > &vec, size_t start, size_t end)
Create an integer interval between start and end.
vector< S > fib
fiber direction
The nodal numbering of the reference mesh (the one stored on HD).
#define SF_MPITAG
the MPI tag when communicating
vector< S > xyz
node cooridnates
T * data()
Pointer to the vector's start.
void read_elements(meshdata< T, S > &mesh, std::string basename)
Read the element data (elements and fibers) of a CARP mesh.
void write_pts_header(FILE *&pts_fd, bool binary, size_t numpts)
Write the header of the points file.
Functions related to network communication.
void write_surface(const meshdata< T, S > &surfmesh, std::string surffile)
void write_fib_block(FILE *&fd, bool write_binary, const meshdata< T, S > &mesh)
Write the local chunk of fibers to a file.
void read_headers(FILE *ele_fd, FILE *fib_fd, bool read_binary, size_t &numelem, bool &twoFib)
Read the header from the element and fiber files.
void insert_points(const vector< S > &pts, const vector< T > &ptsidx, std::list< meshdata< T, S > *> &meshlist)
Insert the points from the read-in buffers into a list of distributed meshes.
void read_pts_block(FILE *&fd, bool read_binary, size_t bsize, vector< S > &xyz)
Read a chunk of points from a file descriptor.
The element numbering of the reference mesh (the one stored on HD).
bool fileExists(std::string filename)
Function which checks if a given file exists.
void read_elem_block(FILE *&fd, bool read_binary, size_t bstart, size_t bsize, meshdata< T, S > &mesh)
Read a block of size bsize from an CARP element file.
size_t read_num_pts(std::string basename)
Function returns the number of points in a CARP points file.
void read_fib_block(FILE *&fd, bool read_binary, bool twoFib, size_t bsize, meshdata< T, S > &mesh)
Read a chunk of fibers from a file descriptor.
elem_t getElemTypeID(char *eletype)
Generate element type enum from string.
void localize(SF_nbr nbr_type)
Localize the connectivity data w.r.t. a given numbering.
hm_int count(const K &key) const
Check if key exists.
size_t size() const
The current size of the vector.
size_t l_numpts
local number of points
A vector storing arbitrary data.
size_t l_numelem
local number of elements
void vec_assign(S *lhs, const V *rhs, size_t size)
Assign the values in rhs to lhs. The data-type of rhs is cast to the type of lhs. ...
vector< S > she
sheet direction
size_t g_numelem
global number of elements
vector< T > & register_numbering(SF_nbr nbr_type)
Register a new numbering to the mesh and return the associated index vector.
void resize(size_t n)
Resize a vector.
vector< elem_t > type
element type
void read_points(const std::string basename, const MPI_Comm comm, vector< S > &pts, vector< T > &ptsidx)
Read the points and insert them into a list of meshes.