1 /*
2 * ***************************************************************************
3 * GAMER = < Geometry-preserving Adaptive MeshER >
4 * Copyright (C) 1994-- Michael Holst and Zeyun Yu
5 *
6 * This library is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU Lesser General Public
8 * License as published by the Free Software Foundation; either
9 * version 2.1 of the License, or (at your option) any later version.
10 *
11 * This library is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with this library; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 * ***************************************************************************
20 */
21
22 /*
23 * ***************************************************************************
24 * File: ReadWriteMesh.C
25 *
26 * Author: Zeyun Yu (zeyun.yu@gmail.com)
27 *
28 * Purpose: Some input and output routines
29 * ***************************************************************************
30 */
31
32 #include <gamer/biom.h>
33 #include <gamer/tetgen.h>
34 #include "gamercf.h"
35
36 /*
37 * All ISO-C headers are coming to us via MALOC, through the includes above.
38 * We now include all architecture-dependent headers.
39 * This buries the architecture-dependent stuff in library, not headers.
40 */
41
42 #if defined(HAVE_UNISTD_H)
43 # include <unistd.h>
44 #endif
45
46 #if defined(HAVE_SYS_TYPES_H)
47 # include <sys/types.h>
48 #endif
49
50 #if defined(HAVE_SYS_TIME_H)
51 # include <sys/time.h>
52 #endif
53
54 #if defined(HAVE_SYS_TIMES_H)
55 # include <sys/times.h>
56 #endif
57
58 #if defined(HAVE_SYS_STAT_H)
59 # include <sys/stat.h>
60 #endif
61
62
63 /*
64 * ***************************************************************************
65 * Routine: same_face
66 *
67 * Author: Johan Hake (hake.dev@gmail.com)
68 *
69 * Purpose: Check if the given numbers relate to the same face
70 * ***************************************************************************
71 */
same_face(int a,int b,int c,int * triface)72 bool same_face(int a, int b, int c, int* triface)
73 {
74 int i;
75 for (i=0; i<3; i++)
76 if (a != triface[i] && b != triface[i] && c != triface[i])
77 return false;
78 return true;
79 }
80
81
82
83 /*
84 * ***************************************************************************
85 * Routine: swap_buffer
86 *
87 * Author: Zeyun Yu (zeyun.yu@gmail.com)
88 *
89 * Purpose: swap the bytes when LITTLE_ENDIAN is enabled
90 * ***************************************************************************
91 */
swap_buffer(char * buffer,int count,int typesize)92 void swap_buffer(char *buffer, int count, int typesize)
93 {
94 char sbuf[4];
95 int i;
96 int temp = 1;
97 unsigned char* chartempf = (unsigned char*) &temp;
98 if(chartempf[0] > '\0') {
99
100 // swapping isn't necessary on single byte data
101 if (typesize == 1)
102 return;
103
104
105 for (i=0; i < count; i++)
106 {
107 memcpy(sbuf, buffer+(i*typesize), typesize);
108
109 switch (typesize)
110 {
111 case 2:
112 {
113 buffer[i*typesize] = sbuf[1];
114 buffer[i*typesize+1] = sbuf[0];
115 break;
116 }
117 case 4:
118 {
119 buffer[i*typesize] = sbuf[3];
120 buffer[i*typesize+1] = sbuf[2];
121 buffer[i*typesize+2] = sbuf[1];
122 buffer[i*typesize+3] = sbuf[0];
123 break;
124 }
125 default:
126 break;
127 }
128 }
129
130 }
131 }
132
133
134
135 /*
136 * ***************************************************************************
137 * Routine: ReadRawiv
138 *
139 * Author: Zeyun Yu (zeyun.yu@gmail.com)
140 *
141 * Purpose: Read the volume in rawiv format
142 (see Chandrajit Bajaj's group for details of this format)
143 * ***************************************************************************
144 */
ReadRawiv(int * xd,int * yd,int * zd,float ** dataset,char * input_name,float * span_t,float * orig_t)145 void ReadRawiv(int *xd, int *yd, int *zd, float **dataset, char *input_name,
146 float *span_t, float *orig_t)
147 {
148 float c_float;
149 unsigned char c_unchar;
150 unsigned short c_unshort;
151 int i,j,k;
152 float *data;
153 int xdim,ydim,zdim;
154 float maxraw;
155 float minraw;
156 float minext[3], maxext[3];
157 int nverts, ncells;
158 unsigned int dim[3];
159 float orig[3], span[3];
160
161 struct stat filestat;
162 size_t size[3];
163 int datatype;
164 int found;
165 FILE *fp;
166 size_t ret;
167
168
169 if ((fp=fopen(input_name, "rb"))==NULL){
170 printf("read error...\n");
171 exit(0);
172 }
173 stat(input_name, &filestat);
174
175 /* reading RAWIV header */
176 ret = fread(minext, sizeof(float), 3, fp);
177 ret = fread(maxext, sizeof(float), 3, fp);
178 ret = fread(&nverts, sizeof(int), 1, fp);
179 ret = fread(&ncells, sizeof(int), 1, fp);
180 #ifdef _LITTLE_ENDIAN
181 swap_buffer((char *)minext, 3, sizeof(float));
182 swap_buffer((char *)maxext, 3, sizeof(float));
183 swap_buffer((char *)&nverts, 1, sizeof(int));
184 swap_buffer((char *)&ncells, 1, sizeof(int));
185 #endif
186
187 size[0] = 12 * sizeof(float) + 2 * sizeof(int) + 3 * sizeof(unsigned int) +
188 nverts * sizeof(unsigned char);
189 size[1] = 12 * sizeof(float) + 2 * sizeof(int) + 3 * sizeof(unsigned int) +
190 nverts * sizeof(unsigned short);
191 size[2] = 12 * sizeof(float) + 2 * sizeof(int) + 3 * sizeof(unsigned int) +
192 nverts * sizeof(float);
193
194 found = 0;
195 for (i = 0; i < 3; i++)
196 if (size[i] == (unsigned int)filestat.st_size)
197 {
198 if (found == 0)
199 {
200 datatype = i;
201 found = 1;
202 }
203 }
204 if (found == 0)
205 {
206 printf("Corrupted file or unsupported dataset type\n");
207 exit(5);
208 }
209
210 ret = fread(dim, sizeof(unsigned int), 3, fp);
211 ret = fread(orig, sizeof(float), 3, fp);
212 ret = fread(span, sizeof(float), 3, fp);
213 #ifdef _LITTLE_ENDIAN
214 swap_buffer((char *)dim, 3, sizeof(unsigned int));
215 swap_buffer((char *)orig, 3, sizeof(float));
216 swap_buffer((char *)span, 3, sizeof(float));
217 #endif
218
219 span_t[0] = span[0];
220 span_t[1] = span[1];
221 span_t[2] = span[2];
222 orig_t[0] = orig[0];
223 orig_t[1] = orig[1];
224 orig_t[2] = orig[2];
225
226 xdim = dim[0];
227 ydim = dim[1];
228 zdim = dim[2];
229
230 data = (float *)malloc(sizeof(float)*xdim*ydim*zdim);
231
232 maxraw = -999999.f;
233 minraw = 999999.f;
234 if (datatype == 0) {
235 printf("data type: unsigned char \n");
236 for (i=0; i<zdim; i++)
237 for (j=0; j<ydim; j++)
238 for (k=0; k<xdim; k++) {
239 ret = fread(&c_unchar, sizeof(unsigned char), 1, fp);
240 data[IndexVect(k,j,i)]=(float)c_unchar;
241 if (c_unchar > maxraw)
242 maxraw = c_unchar;
243 if (c_unchar < minraw)
244 minraw = c_unchar;
245 }
246 }
247 else if (datatype == 1) {
248 printf("data type: unsigned short \n");
249 for (i=0; i<zdim; i++)
250 for (j=0; j<ydim; j++)
251 for (k=0; k<xdim; k++) {
252 ret = fread(&c_unshort, sizeof(unsigned short), 1, fp);
253 #ifdef _LITTLE_ENDIAN
254 swap_buffer((char *)&c_unshort, 1, sizeof(unsigned short));
255 #endif
256 data[IndexVect(k,j,i)]=(float)c_unshort;
257
258 if (c_unshort > maxraw)
259 maxraw = c_unshort;
260 if (c_unshort < minraw)
261 minraw = c_unshort;
262 }
263 }
264 else if (datatype == 2) {
265 printf("data type: float \n");
266 for (i=0; i<zdim; i++)
267 for (j=0; j<ydim; j++)
268 for (k=0; k<xdim; k++) {
269 ret = fread(&c_float, sizeof(float), 1, fp);
270 #ifdef _LITTLE_ENDIAN
271 swap_buffer((char *)&c_float, 1, sizeof(float));
272 #endif
273 data[IndexVect(k,j,i)]=c_float;
274
275 if (c_float > maxraw)
276 maxraw = c_float;
277 if (c_float < minraw)
278 minraw = c_float;
279 }
280 }
281 else {
282 printf("error\n");
283 fclose(fp);
284 exit(1);
285 }
286
287 fclose(fp);
288 printf("minimum = %f, maximum = %f \n",minraw,maxraw);
289
290 for (i=0; i<zdim; i++)
291 for (j=0; j<ydim; j++)
292 for (k=0; k<xdim; k++) {
293 data[IndexVect(k,j,i)] = 255*(data[IndexVect(k,j,i)] -
294 minraw)/(maxraw-minraw);
295 }
296
297 printf("dimension: %d X %d X %d\n",xdim,ydim,zdim);
298 *xd = xdim;
299 *yd = ydim;
300 *zd = zdim;
301 *dataset = data;
302
303 }
304
305
306
307 /*
308 * ***************************************************************************
309 * SubRoutine: rgb_to_marker
310 *
311 * Author: Johan Hake (hake.dev@gmail.com)
312 *
313 * Purpose: Transform an rgb value to an integer
314 *
315 * ***************************************************************************
316 */
rgb_to_marker(float r,float g,float b)317 int rgb_to_marker(float r, float g, float b)
318 {
319 if (r < 0 | r > 1 | g < 0 | g > 1 | b < 0 | b > 1 ){
320 printf("Expected individual RGB value to be betwen 0 and 1.\n");
321 exit(1);
322 }
323
324 return round(r*10)*121 + round(g*10)*11 + round(b*10);
325 }
326
327
328 /*
329 * ***************************************************************************
330 * Routine: SurfaceMesh_readOFF
331 *
332 * Author: Zeyun Yu (zeyun.yu@gmail.com)
333 *
334 * Purpose: Read a surface mesh in OFF format, if a volumetric mesh is
335 * provided, a surface mesh will be extracted
336 * ***************************************************************************
337 */
SurfaceMesh_readOFF(char * input_name)338 SurfaceMesh* SurfaceMesh_readOFF(char *input_name)
339 {
340 unsigned int n,m;
341 unsigned int a,b,c,d;
342 float x,y,z, color_r, color_g, color_b, color_a;
343 TeTraMesh *volmesh = NULL;
344 int v_n, t_n, e_n, character;
345 char line[256];
346 FILE *fin = NULL;
347 fpos_t position;
348 int* face_markers = NULL;
349 bool has_marker = false;
350 SurfaceMesh* surfmesh = NULL;
351
352
353 if ((fin = fopen(input_name, "r"))==NULL){
354 printf("Read error. File \'%s\' could not be read.\n", input_name);
355 exit(1);
356 };
357
358 while (fgets(line, 256, fin) != NULL) {
359 if (line[0]=='O' && line[1]=='F' && line[2]=='F')
360 break;
361 }
362
363 if (fscanf(fin,"%d %d %d\n", &v_n, &t_n, &e_n) != 3){
364 printf("Read error. Expected 3 integer number for the number of vertices "
365 "and simplices.\n");
366 exit(1);
367 }
368
369 printf(" vertices: %d --- simplices: %d \n", v_n, t_n);
370
371 // Assign memory
372 surfmesh = SurfaceMesh_ctor(v_n, t_n);
373
374 // Read vertex coordinates
375 for (n = 0; n < surfmesh->nv; n++) {
376 if (fscanf(fin, "%f %f %f\n", &x, &y, &z) != 3){
377 printf("Read error. Expected 3 floats for the coordinates of vertex: %d.\n", n);
378 exit(1);
379 }
380
381 surfmesh->vertex[n].x = x;
382 surfmesh->vertex[n].y = y;
383 surfmesh->vertex[n].z = z;
384 }
385
386 // Read the number of vertices per simplex from the first line of simplices
387 n = fscanf(fin, "%d", &m);
388
389 // Check format of the read simplices
390 if (m != 3 && m != 4) {
391 printf("Read error. Expected a 3 or 4 for the first value in the first simplex line.\n");
392 exit(1);
393 }
394
395 // Input is surface mesh
396 else if (m == 3) {
397 printf(" Input is surface mesh.\n");
398
399 // Read the rest of the first line
400 if (fscanf(fin, "%d %d %d", &a, &b, &c) != 3){
401 printf("Read error. Expected 3 integers for the first simplex.\n", n);
402 exit(1);
403 }
404
405 // Get position of file and check if we have more symbols
406 fgetpos(fin, &position);
407
408 // Read any blanks
409 character = fgetc(fin);
410 while (character == ' '){
411 character = fgetc(fin);
412 }
413
414 // If we do not have en end of line character we have marker
415 if (character != '\n'){
416
417 // Set marker flag
418 has_marker = true;
419
420 // Rewind the file pointer
421 fsetpos(fin, &position);
422
423 // Reset facemarkers
424 SurfaceMesh_resetFaceMarkers(surfmesh);
425
426 // Read first rgb values
427 if (fscanf(fin, "%f %f %f %f", &color_r, &color_g, &color_b, &color_a) != 4){
428 printf("Read error. Expected 4 floats for the RGBA values of the first face.\n");
429 exit(1);
430 }
431
432 // Get the first face marker
433 surfmesh->face_markers[0] = rgb_to_marker(color_r, color_g, color_b);
434
435 while (fgetc(fin) != '\n'){}
436 }
437
438 // Assign first simplex
439 surfmesh->face[0].a = a;
440 surfmesh->face[0].b = b;
441 surfmesh->face[0].c = c;
442
443 // Read the rest of the simplices
444 for (n = 1; n < surfmesh->nf; n++) {
445 if (fscanf(fin, "%d %d %d %d", &m, &a, &b, &c) != 4){
446 printf("Read error. Expected 4 integers for simplex %d.\n", n);
447 exit(1);
448 }
449 surfmesh->face[n].a = a;
450 surfmesh->face[n].b = b;
451 surfmesh->face[n].c = c;
452
453 // If we have face markers
454 if (has_marker){
455
456 // Read first rgb values
457 if (fscanf(fin, "%f %f %f %f", &color_r, &color_g, &color_b, &color_a) != 4){
458 printf("Read error. Expected 4 floats for the RGBA values of face: %d.\n", n);
459 exit(1);
460 }
461
462 // Get the other face markers
463 surfmesh->face_markers[n] = rgb_to_marker(color_r, color_g, color_b);
464
465 }
466
467 // Skip any additional character on this line
468 while (fgetc(fin) != '\n'){}
469 }
470
471 // Close file
472 fclose(fin);
473 }
474
475 // Input is tetrahedral mesh.
476 else {
477 printf(" Input is volumetric mesh...\n");
478 fclose(fin);
479
480 // Create tetrahedal mesh
481 volmesh = (TeTraMesh*)malloc(sizeof(TeTraMesh));
482 volmesh->nv = surfmesh->nv;
483 volmesh->nf = surfmesh->nf;
484 volmesh->neighbor = NULL;
485
486 // We have allready read the vertex coordinates
487 volmesh->vertex = surfmesh->vertex;
488
489 // Assign memory for faces
490 volmesh->face = (INT4VECT *)malloc(sizeof(INT4VECT)*volmesh->nf);
491
492 // Read the rest of the first line
493 if (fscanf(fin, "%d %d %d %d", &a, &b, &c, &d) != 4){
494 printf("Read error. Expected 4 integers for the first simplex.\n", n);
495 exit(1);
496 }
497
498 // Skip any additional character on this line
499 while (fgetc(fin) != '\n'){}
500
501 // Assign first simplex from allready read data
502 volmesh->face[0].a = a;
503 volmesh->face[0].b = b;
504 volmesh->face[0].c = c;
505 volmesh->face[0].d = d;
506
507 // Read the rest of the simplices
508 for (n = 1; n < volmesh->nf; n++) {
509 if (fscanf(fin, "%d %d %d %d %d", &m, &a, &b, &c, &d) != 5){
510 printf("Read error. Expected 5 integers for simplex %d.\n", n);
511 exit(1);
512 }
513 volmesh->face[n].a = a;
514 volmesh->face[n].b = b;
515 volmesh->face[n].c = c;
516 volmesh->face[n].d = d;
517
518 // Skip any additional character on this line
519 while (fgetc(fin) != '\n'){}
520 }
521
522 fclose(fin);
523
524 SurfaceExtract(volmesh, surfmesh);
525 free(volmesh->vertex);
526 free(volmesh->face);
527 free(volmesh);
528
529 }
530
531 return surfmesh;
532 }
533
534
535
536
537 /*
538 * ***************************************************************************
539 * Routine: write_rawiv_float
540 *
541 * Author: Zeyun Yu (zeyun.yu@gmail.com)
542 *
543 * Purpose: write the volume into a rawiv file
544 * ***************************************************************************
545 */
write_rawiv_float(FILE * fp,float * dataset,int dim[3],float min[3],float max[3])546 void write_rawiv_float(FILE* fp,float *dataset,int dim[3],float min[3],float max[3])
547 {
548 int i, j, k;
549 float c_float;
550 float orig[3],span[3];
551 int xdim,ydim,zdim;
552 int nverts,ncells;
553
554
555 xdim = dim[0];
556 ydim = dim[1];
557 zdim = dim[2];
558
559 orig[0] = min[0];
560 orig[1] = min[1];
561 orig[2] = min[2];
562 span[0] = (max[0] - min[0]) / (float)(xdim-1);
563 span[1] = (max[1] - min[1]) / (float)(ydim-1);
564 span[2] = (max[2] - min[2]) / (float)(zdim-1);
565 nverts = xdim*ydim*zdim;
566 ncells = (xdim-1)*(ydim-1)*(zdim-1);
567
568
569 #ifdef _LITTLE_ENDIAN
570 swap_buffer((char *)min, 3, sizeof(float));
571 swap_buffer((char *)max, 3, sizeof(float));
572 swap_buffer((char *)&nverts, 1, sizeof(int));
573 swap_buffer((char *)&ncells, 1, sizeof(int));
574 swap_buffer((char *)dim, 3, sizeof(unsigned int));
575 swap_buffer((char *)orig, 3, sizeof(float));
576 swap_buffer((char *)span, 3, sizeof(float));
577 #endif
578
579 fwrite(min, sizeof(float), 3, fp);
580 fwrite(max, sizeof(float), 3, fp);
581 fwrite(&nverts, sizeof(int), 1, fp);
582 fwrite(&ncells, sizeof(int), 1, fp);
583 fwrite(dim, sizeof(unsigned int), 3, fp);
584 fwrite(orig, sizeof(float), 3, fp);
585 fwrite(span, sizeof(float), 3, fp);
586 #ifdef _LITTLE_ENDIAN
587 swap_buffer((char *)min, 3, sizeof(float));
588 swap_buffer((char *)max, 3, sizeof(float));
589 swap_buffer((char *)&nverts, 1, sizeof(int));
590 swap_buffer((char *)&ncells, 1, sizeof(int));
591 swap_buffer((char *)dim, 3, sizeof(unsigned int));
592 swap_buffer((char *)orig, 3, sizeof(float));
593 swap_buffer((char *)span, 3, sizeof(float));
594 #endif
595
596 for (k=0; k<zdim; k++)
597 for (j=0; j<ydim; j++)
598 for (i=0; i<xdim; i++) {
599 c_float = dataset[IndexVect(i,j,k)];
600 #ifdef _LITTLE_ENDIAN
601 swap_buffer((char *)&c_float, 1, sizeof(float));
602 #endif
603 fwrite(&c_float, sizeof(float), 1, fp);
604 #ifdef _LITTLE_ENDIAN
605 swap_buffer((char *)&c_float, 1, sizeof(float));
606 #endif
607 }
608
609 fclose(fp);
610 }
611
612
613 /*
614 * ***************************************************************************
615 * Routine: WriteSurfaceMeshOFF
616 *
617 * Author: Zeyun Yu (zeyun.yu@gmail.com)
618 *
619 * Purpose: write a surface mesh into an OFF file
620 * ***************************************************************************
621 */
SurfaceMesh_writeOFF(SurfaceMesh * surfmesh,char * filename)622 void SurfaceMesh_writeOFF(SurfaceMesh *surfmesh, char* filename)
623 {
624 int i;
625 FILE* fout;
626
627 // Try to open the file
628 if ((fout=fopen(filename, "wb"))==NULL){
629 printf("write error...\n");
630 exit(0);
631 };
632
633 fprintf(fout, "OFF\n");
634 fprintf(fout, "%d %d %d\n",surfmesh->nv,surfmesh->nf,surfmesh->nv+surfmesh->nf-2);
635
636 for (i = 0; i < surfmesh->nv; i++)
637 fprintf(fout, "%17.10e %17.10e %17.10e\n",surfmesh->vertex[i].x,surfmesh->vertex[i].y,surfmesh->vertex[i].z);
638
639 for (i = 0; i < surfmesh->nf; i++) {
640 fprintf(fout, "3 %d %d %d\n",surfmesh->face[i].a,surfmesh->face[i].b,surfmesh->face[i].c);
641 }
642 fclose(fout);
643 }
644
645 /*
646 * ***************************************************************************
647 * Routine: write_off
648 *
649 * Author: Zeyun Yu (zeyun.yu@gmail.com)
650 *
651 * Purpose: Write a surface mesh into an OFF format
652 * ***************************************************************************
653 */
GemMesh_writeOFF(GemMesh * Gem_mesh,char * filename)654 void GemMesh_writeOFF(GemMesh* Gem_mesh, char* filename)
655 {
656 long i;
657
658 FILE* fout;
659
660 // Try to open the file
661 if ((fout=fopen(filename, "wb"))==NULL){
662 printf("write error...\n");
663 exit(0);
664 };
665
666 fprintf(fout, "OFF\n");
667 fprintf(fout, "%d %d 0\n", Gem_mesh->vertices, Gem_mesh->simplices);
668 for (i = 0; i < Gem_mesh->vertices; i++)
669 fprintf(fout, "%.12f %.12f %.12f\n", Gem_mesh->vv[i].x, Gem_mesh->vv[i].y, Gem_mesh->vv[i].z);
670
671 for (i = 0; i < Gem_mesh->simplices; i++)
672 fprintf(fout, "4 %8d %8d %8d %8d\n", Gem_mesh->ss[i].na, Gem_mesh->ss[i].nb,
673 Gem_mesh->ss[i].nc, Gem_mesh->ss[i].nd);
674
675 fclose(fout);
676 }
677
678 /*
679 * ***************************************************************************
680 * Routine: GemMesh_dtor
681 *
682 * Author: Johan Hake (hake.dev@gmail.com)
683 *
684 * Purpose: Release memory from a GemMesh
685 * ***************************************************************************
686 */
687 // FIXME: Move all GemMesh functions to its own
GemMesh_dtor(GemMesh * gem_mesh)688 void GemMesh_dtor(GemMesh* gem_mesh)
689 {
690 free(gem_mesh->vv);
691 free(gem_mesh->ss);
692 free(gem_mesh);
693 }
694
695 /*
696 * ***************************************************************************
697 * Routine: GemMesh_fromSurfaceMesh
698 *
699 * Author: Johan Hake (hake.dev@gmail.com)
700 *
701 * Purpose: Use tetgen to generate a GemMesh from a surface mesh
702 * ***************************************************************************
703 */
GemMesh_fromSurfaceMesh(SurfaceMesh * surfmesh,char * tetgen_params)704 GemMesh* GemMesh_fromSurfaceMesh(SurfaceMesh* surfmesh, char* tetgen_params)
705 {
706
707 unsigned int j;
708
709 tetgenio in, out;
710 tetgenio::facet *f;
711 tetgenio::polygon *p;
712 bool use_facemarkers = false;
713 in.firstnumber = 1;
714
715 // Assign vertex information
716 in.numberofpoints = surfmesh->nv;
717 in.pointlist = new REAL[in.numberofpoints * 3];
718 for (j = 0; j < in.numberofpoints; j++) {
719 in.pointlist[j*3+0] = surfmesh->vertex[j].x;
720 in.pointlist[j*3+1] = surfmesh->vertex[j].y;
721 in.pointlist[j*3+2] = surfmesh->vertex[j].z;
722 }
723
724 // Assign face information
725 use_facemarkers = surfmesh->nfm == surfmesh->nf;
726 in.numberoffacets = surfmesh->nf;
727 in.facetlist = new tetgenio::facet[in.numberoffacets];
728 if (use_facemarkers)
729 in.facetmarkerlist = new int[surfmesh->nf];
730
731 for (j = 0; j < in.numberoffacets; j++) {
732 f = &in.facetlist[j];
733 f->holelist = (REAL *) NULL;
734 f->numberofholes = 0;
735 f->numberofpolygons = 1;
736 f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
737 p = &f->polygonlist[0];
738 p->numberofvertices = 3;
739 p->vertexlist = new int[p->numberofvertices];
740 p->vertexlist[0] = surfmesh->face[j].a+in.firstnumber;
741 p->vertexlist[1] = surfmesh->face[j].b+in.firstnumber;
742 p->vertexlist[2] = surfmesh->face[j].c+in.firstnumber;
743 if (use_facemarkers)
744 in.facetmarkerlist[j] = surfmesh->face_markers[j];
745 }
746
747 // Add boundary marker on each node
748 // FIXME: Why? Aren't the markers set on the generatedmesh?
749 in.pointmarkerlist = new int[in.numberofpoints];
750 for (j = 0; j < in.numberofpoints; j++)
751 in.pointmarkerlist[j] = 1;
752
753 // Debug
754 in.save_nodes("plc");
755 in.save_poly("plc");
756
757 // Call TetGen
758 tetrahedralize(tetgen_params, &in, &out, NULL);
759
760 // Debug
761 out.save_nodes("result");
762 out.save_elements("result");
763 out.save_faces("result");
764
765 // Convert to a generalized tetrahedral mesh structure
766 return GemMesh_fromTetgen(out);
767
768 }
769
770
771 /*
772 * ***************************************************************************
773 * Routine: GemMesh_fromTetgen
774 *
775 * Author: Johan Hake (hake.dev@gmail.com)
776 *
777 * Purpose: Generate a GemMesh structure
778 * ***************************************************************************
779 */
780 /*
781 FIXME: Are we using the neighbor list at all?
782 */
GemMesh_fromTetgen(tetgenio & tetio)783 GemMesh* GemMesh_fromTetgen(tetgenio& tetio)
784 {
785 GemMesh *Gem_mesh;
786 unsigned int i, j, k;
787 FILE *fout;
788 long node_index;
789 int *face_type, tetra_node[4];
790 float x,y,z;
791 int cell_type;
792 float dist;
793 long a,b,c,d;
794 float ax,ay,az;
795 float bx,by,bz;
796 float cx,cy,cz;
797 float dx,dy,dz;
798 int active, exterior_cell_type;
799 long neighbor, neightype;
800 long *map_w_o,*map_w_i,*map_w_t;
801 Triface** tet_to_triface;
802 int tet0, tet1, tet2;
803 int* vertex_markers;
804 Triface* triface;
805
806 // Initialization and memory allocation
807 Gem_mesh = (GemMesh*)malloc(sizeof(GemMesh));
808 Gem_mesh->dim = 3;
809 Gem_mesh->dimii = 3;
810 Gem_mesh->vv = (FETK_VX *)malloc(sizeof(FETK_VX)*tetio.numberofpoints);
811 Gem_mesh->ss = (FETK_SS *)malloc(sizeof(FETK_SS)*tetio.numberoftetrahedra);
812
813 tet_to_triface = new Triface*[tetio.numberoftetrahedra];
814 vertex_markers = new int[tetio.numberofpoints];
815
816 for (i = 0; i < tetio.numberoftetrahedra; i++){
817 tet_to_triface[i] = NULL;
818 }
819
820 for (i = 0; i < tetio.numberofpoints; i++){
821 vertex_markers[i] = 0;
822 }
823
824 // Using the markers from the trifaces to mark vertices beeing on the boundary
825 for (i = 0; i < tetio.numberoftrifaces; i++){
826 tet0 = tetio.adjtetlist[i*2]-1;
827 tet1 = tetio.adjtetlist[i*2+1]-1;
828 tet2 = (tet0 == -2) ? tet1 : ((tet1 == -2) ? tet0 : -2);
829 // Get the tet to triface info
830 if (tet2 != -2){
831 triface = new Triface;
832 triface->next = NULL;
833 triface->points[0] = tetio.trifacelist[i*3]-1;
834 triface->points[1] = tetio.trifacelist[i*3+1]-1;
835 triface->points[2] = tetio.trifacelist[i*3+2]-1;
836 if (tetio.trifacemarkerlist != NULL)
837 triface->marker = tetio.trifacemarkerlist[i];
838 else
839 triface->marker = 1;
840 if (tet_to_triface[tet2] != NULL)
841 triface->next = tet_to_triface[tet2];
842 tet_to_triface[tet2] = triface;
843 vertex_markers[triface->points[0]] = triface->marker;
844 vertex_markers[triface->points[1]] = triface->marker;
845 vertex_markers[triface->points[2]] = triface->marker;
846 }
847 }
848
849 // Generate facetype information
850 face_type = new int[tetio.numberoftetrahedra*4];
851
852 for (i = 0; i < tetio.numberoftetrahedra; i++) {
853
854 // Zero out all face types
855 for (j = 0; j < tetio.numberofcorners; j++)
856 face_type[i*4+j] = 0;
857
858 // Check tetra orientation
859 a = tetio.tetrahedronlist[i * tetio.numberofcorners]-1;
860 b = tetio.tetrahedronlist[i * tetio.numberofcorners + 1]-1;
861 c = tetio.tetrahedronlist[i * tetio.numberofcorners + 2]-1;
862 d = tetio.tetrahedronlist[i * tetio.numberofcorners + 3]-1;
863 ax = tetio.pointlist[a * 3];
864 ay = tetio.pointlist[a * 3 + 1];
865 az = tetio.pointlist[a * 3 + 2];
866 bx = tetio.pointlist[b * 3];
867 by = tetio.pointlist[b * 3 + 1];
868 bz = tetio.pointlist[b * 3 + 2];
869 cx = tetio.pointlist[c * 3];
870 cy = tetio.pointlist[c * 3 + 1];
871 cz = tetio.pointlist[c * 3 + 2];
872 dx = tetio.pointlist[d * 3];
873 dy = tetio.pointlist[d * 3 + 1];
874 dz = tetio.pointlist[d * 3 + 2];
875 bx -= ax;
876 by -= ay;
877 bz -= az;
878 cx -= ax;
879 cy -= ay;
880 cz -= az;
881 dx -= ax;
882 dy -= ay;
883 dz -= az;
884 ax = by*cz-bz*cy;
885 ay = bz*cx-bx*cz;
886 az = bx*cy-by*cx;
887 if (ax*dx+ay*dy+az*dz < 0) {
888 tetio.tetrahedronlist[i * tetio.numberofcorners + 1] = c;
889 tetio.tetrahedronlist[i * tetio.numberofcorners + 2] = b;
890 k = tetio.neighborlist[i * 4 + 1];
891 tetio.neighborlist[i * 4 + 1] = tetio.neighborlist[i * 4 + 2];
892 tetio.neighborlist[i * 4 + 2] = k;
893 }
894
895 // Check if we are on the boundary
896 while (tet_to_triface[i] != NULL){
897
898 // If so find out which face is on the boundary
899 if (same_face(b, c, d, tet_to_triface[i]->points)){
900 face_type[i*4] = tet_to_triface[i]->marker;
901 triface = tet_to_triface[i];
902 tet_to_triface[i] = triface->next;
903 delete triface;
904 continue;
905 }
906
907 if (same_face(a, c, d, tet_to_triface[i]->points)){
908 face_type[i*4+1] = tet_to_triface[i]->marker;
909 triface = tet_to_triface[i];
910 tet_to_triface[i] = triface->next;
911 delete triface;
912 continue;
913 }
914
915 if (same_face(a, b, d, tet_to_triface[i]->points)){
916 face_type[i*4+2] = tet_to_triface[i]->marker;
917 triface = tet_to_triface[i];
918 tet_to_triface[i] = triface->next;
919 delete triface;
920 continue;
921 }
922
923 if (same_face(a, b, c, tet_to_triface[i]->points)){
924 face_type[i*4+3] = tet_to_triface[i]->marker;
925 triface = tet_to_triface[i];
926 tet_to_triface[i] = triface->next;
927 delete triface;
928 continue;
929 }
930 }
931 }
932
933 // Start ouput mesh
934 Gem_mesh->vertices = tetio.numberofpoints;
935 Gem_mesh->simplices = tetio.numberoftetrahedra;
936
937 // Write all nodes to a GemMesh structure
938 for (i = 0; i < tetio.numberofpoints; i++){
939 Gem_mesh->vv[i].id = i;
940 Gem_mesh->vv[i].chrt = vertex_markers[i]; // This might be ambigious
941 Gem_mesh->vv[i].x = tetio.pointlist[i * 3];
942 Gem_mesh->vv[i].y = tetio.pointlist[i * 3 + 1];
943 Gem_mesh->vv[i].z = tetio.pointlist[i * 3 + 2];
944 }
945
946 // Write all tets to a GemMesh structure
947 for (i = 0; i < tetio.numberoftetrahedra; i++) {
948 Gem_mesh->ss[i].id = i;
949 Gem_mesh->ss[i].grp = 0;
950 if (tetio.numberoftetrahedronattributes>0)
951 Gem_mesh->ss[i].mat = (int)tetio.tetrahedronattributelist[i*tetio.numberoftetrahedronattributes];
952 else
953 Gem_mesh->ss[i].mat = 0;
954 Gem_mesh->ss[i].fa = face_type[4*i];
955 Gem_mesh->ss[i].fb = face_type[4*i+1];
956 Gem_mesh->ss[i].fc = face_type[4*i+2];
957 Gem_mesh->ss[i].fd = face_type[4*i+3];
958
959 // NOTE: Index in tetrahedronlist starts from 1
960 Gem_mesh->ss[i].na = tetio.tetrahedronlist[i*tetio.numberofcorners] - 1;
961 Gem_mesh->ss[i].nb = tetio.tetrahedronlist[i*tetio.numberofcorners+1] - 1;
962 Gem_mesh->ss[i].nc = tetio.tetrahedronlist[i*tetio.numberofcorners+2] - 1;
963 Gem_mesh->ss[i].nd = tetio.tetrahedronlist[i*tetio.numberofcorners+3] - 1;
964 }
965
966 // Release memory
967 delete[] face_type;
968 delete[] tet_to_triface;
969 delete[] vertex_markers;
970
971 return Gem_mesh;
972 }
973
974 /*
975 * ***************************************************************************
976 * Routine: GemMesh_writeMcsf
977 *
978 * Author: Johan Hake (hake.dev@gmail.com)
979 *
980 * Purpose: Write tetrahedralized mesh to disk (MCSF formats)
981 * ***************************************************************************
982 */
983
GemMesh_writeMcsf(GemMesh * out,char * filename)984 void GemMesh_writeMcsf(GemMesh* out, char* filename)
985 {
986 unsigned int i;
987 FILE *fout;
988 FETK_VX* vv;
989 FETK_SS* ss;
990
991 // Open file
992 if ((fout=fopen(filename, "wb"))==NULL){
993 printf("write error...\n");
994 exit(0);
995 };
996
997 // File head
998 fprintf(fout,"mcsf_begin=1; \n \n");
999 fprintf(fout," dim = %d; \n", out->dim);
1000 fprintf(fout," dimii = %d; \n", out->dimii);
1001 fprintf(fout," vertices = %d;\n", out->vertices);
1002 fprintf(fout,"simplices = %d;\n", out->simplices);
1003
1004 // Node head
1005 fprintf(fout,"vert=[ \n");
1006 fprintf(fout,"%%------------------------------------------------------------------------------ \n");
1007 fprintf(fout,"%% Node-ID Chrt X-Coordinate Y-coordinate Z-coordinate \n");
1008 fprintf(fout,"%%--------- ---- ---------------- ---------------- ---------------- \n");
1009
1010 // Write all nodes with type to file
1011 for (i = 0; i < out->vertices; i++){
1012 vv = out->vv + i;
1013 fprintf(fout, "%10d %d %17.10e %17.10e %17.10e \n", i,
1014 vv->chrt, vv->x, vv->y, vv->z);
1015 }
1016
1017 fprintf(fout,"]; \n");
1018
1019 // Cell (simplex) head
1020 fprintf(fout,"simp=[ \n");
1021 fprintf(fout,"%%------------------------------------------------------------------------------------------ \n");
1022 fprintf(fout,"%% Simp-ID Grp Mat Face-Types Vertex-Numbers \n");
1023 fprintf(fout,"%%--------- --- --- --------------------- ------------------------------------------- \n");
1024
1025 for (i = 0; i < out->simplices; i++) {
1026 ss = out->ss + i;
1027 fprintf(fout,"%10d 0 %5d %5d %5d %5d %5d %10d %10d %10d %10d \n", i, ss->id,
1028 ss->fa, ss->fb, ss->fc, ss->fd, ss->na, ss->nb, ss->nc, ss->nd);
1029 }
1030
1031 fprintf(fout,"]; \n");
1032 fprintf(fout,"mcsf_end=1; \n");
1033 fclose(fout);
1034
1035 }
1036
1037 /*
1038 * ***************************************************************************
1039 * Routine: write_mcsf
1040 *
1041 * Author: Zeyun Yu (zeyun.yu@gmail.com)
1042 *
1043 * Purpose: Write surface/tetrahedral meshes into disks (MCSF/OFF formats)
1044 * ***************************************************************************
1045 */
1046 //void tetgenio::write_mcsf(char* filename, unsigned int* ActiveSite)
1047 //{
1048 // unsigned int i, j, k;
1049 // FILE *fout;
1050 // int *node_attr;
1051 // long node_index;
1052 // int *face_type, tetra_node[4];
1053 // float x,y,z;
1054 // char file_name[256];
1055 // int cell_type, interior_node_type, exterior_node_type;
1056 // float dist;
1057 // long a,b,c,d;
1058 // float ax,ay,az;
1059 // float bx,by,bz;
1060 // float cx,cy,cz;
1061 // float dx,dy,dz;
1062 // int active, exterior_cell_type;
1063 // long neighbor, neightype;
1064 // long *map_w_o,*map_w_i,*map_w_t;
1065 //
1066 // printf("Num trifaces: %d\n", numberoftrifaces);
1067 //
1068 // // Check the type of exterior tetrahedra
1069 // // In Tetgen, it could be radmonly 1 or 2 (or even other integers if more materials are present)
1070 // fflush(stdout);
1071 // exterior_cell_type = -1;
1072 // for (i = 0; i < numberoftetrahedra; i++) {
1073 // cell_type = (int)tetrahedronattributelist[i * numberoftetrahedronattributes];
1074 // for (j = 0; j < numberofcorners; j++) {
1075 // neighbor = neighborlist[i * 4 + j] - 1;
1076 // if (neighbor < 0) {
1077 // exterior_cell_type = cell_type;
1078 // break;
1079 // }
1080 // }
1081 // if (exterior_cell_type >= 0)
1082 // break;
1083 // }
1084 // printf("exterior_cell_type = %d \n", exterior_cell_type);
1085 //
1086 // // assign node attributes
1087 // // Node attributes are assigned as either being situated at the exterior or
1088 // // interior domain
1089 // // For now we go for:
1090 // // 0 = interior domain
1091 // // 1 = exterior domain
1092 // // This might be user specific in future
1093 // interior_node_type = 0;
1094 // exterior_node_type = 1;
1095 //
1096 // node_attr = new int[numberofpoints];
1097 // for (i = 0; i < numberofpoints; i++) {
1098 // node_attr[i] = 0;
1099 // }
1100 //
1101 // for (i = 0; i < numberoftetrahedra; i++) {
1102 // cell_type = (int)tetrahedronattributelist[i * numberoftetrahedronattributes];
1103 // for (j = 0; j < numberofcorners; j++) {
1104 // // Note that index in tetrahedronlist starts from 1 !!!!!!!!
1105 // node_index = tetrahedronlist[i * numberofcorners + j] - 1;
1106 //
1107 // // If exterior_cell_type is 1
1108 // if (exterior_cell_type == 1) {
1109 // if (cell_type >= 3) cell_type = 2;
1110 // // FIXME: I do not understand this logic. Slowly getting there... JH
1111 // node_attr[node_index] = cell_type==1 ? exterior_node_type : interior_node_type;
1112 // }
1113 // else if (exterior_cell_type == 2) {
1114 // if (cell_type >= 3) cell_type = 1;
1115 // node_attr[node_index] = cell_type==2 ? exterior_node_type : interior_node_type;
1116 // }
1117 // else {
1118 // printf("there are extra materials inside molecular surfaces\n");
1119 // exit(0);
1120 // }
1121 // }
1122 //
1123 // // check tetra orientation
1124 // a = tetrahedronlist[i * numberofcorners]-1;
1125 // b = tetrahedronlist[i * numberofcorners + 1]-1;
1126 // c = tetrahedronlist[i * numberofcorners + 2]-1;
1127 // d = tetrahedronlist[i * numberofcorners + 3]-1;
1128 // ax = pointlist[a * 3];
1129 // ay = pointlist[a * 3 + 1];
1130 // az = pointlist[a * 3 + 2];
1131 // bx = pointlist[b * 3];
1132 // by = pointlist[b * 3 + 1];
1133 // bz = pointlist[b * 3 + 2];
1134 // cx = pointlist[c * 3];
1135 // cy = pointlist[c * 3 + 1];
1136 // cz = pointlist[c * 3 + 2];
1137 // dx = pointlist[d * 3];
1138 // dy = pointlist[d * 3 + 1];
1139 // dz = pointlist[d * 3 + 2];
1140 // bx -= ax;
1141 // by -= ay;
1142 // bz -= az;
1143 // cx -= ax;
1144 // cy -= ay;
1145 // cz -= az;
1146 // dx -= ax;
1147 // dy -= ay;
1148 // dz -= az;
1149 // ax = by*cz-bz*cy;
1150 // ay = bz*cx-bx*cz;
1151 // az = bx*cy-by*cx;
1152 // if (ax*dx+ay*dy+az*dz < 0) {
1153 // tetrahedronlist[i * numberofcorners + 1] = c;
1154 // tetrahedronlist[i * numberofcorners + 2] = b;
1155 // k = neighborlist[i * 4 + 1];
1156 // neighborlist[i * 4 + 1] = neighborlist[i * 4 + 2];
1157 // neighborlist[i * 4 + 2] = k;
1158 // }
1159 // }
1160 //
1161 // // output the whole tetra mesh
1162 // if ((fout=fopen(filename, "wb"))==NULL){
1163 // printf("write error...\n");
1164 // exit(0);
1165 // };
1166 //
1167 // // File head
1168 // fprintf(fout,"mcsf_begin=1; \n \n");
1169 // fprintf(fout," dim = 3; \n");
1170 // fprintf(fout," dimii = 3; \n");
1171 // fprintf(fout," vertices = %d;\n", numberofpoints);
1172 // fprintf(fout,"simplices = %d;\n", numberoftetrahedra);
1173 //
1174 // // Node head
1175 // fprintf(fout,"vert=[ \n");
1176 // fprintf(fout,"%%------------------------------------------------------------------------------ \n");
1177 // fprintf(fout,"%% Node-ID Chrt X-Coordinate Y-coordinate Z-coordinate \n");
1178 // fprintf(fout,"%%--------- ---- ---------------- ---------------- ---------------- \n");
1179 //
1180 // // Write all nodes with type to file
1181 // for (i = 0; i < numberofpoints; i++)
1182 // fprintf(fout, "%10d %d %17.10e %17.10e %17.10e \n", i,
1183 // node_attr[i], pointlist[i*3], pointlist[i*3+1], pointlist[i*3+2]);
1184 //
1185 // fprintf(fout,"]; \n");
1186 //
1187 // // Cell (simplex) head
1188 // fprintf(fout,"simp=[ \n");
1189 // fprintf(fout,"%%------------------------------------------------------------------------------------------ \n");
1190 // fprintf(fout,"%% Simp-ID Grp Mat Face-Types Vertex-Numbers \n");
1191 // fprintf(fout,"%%--------- --- --- --------------------- ------------------------------------------- \n");
1192 //
1193 // face_type = new int[numberoftetrahedra*4];
1194 // for (i = 0; i < numberoftetrahedra; i++) {
1195 // for (j = 0; j < numberofcorners; j++)
1196 // // NOTE: Index in tetrahedronlist starts from 1
1197 // tetra_node[j] = tetrahedronlist[i * numberofcorners + j] - 1;
1198 //
1199 // cell_type = (int)tetrahedronattributelist[i * numberoftetrahedronattributes];
1200 //
1201 // // Check face type
1202 // for (j = 0; j < numberofcorners; j++) {
1203 // face_type[4*i+j] = 0;
1204 //
1205 // // NOTE: Index in neighborlist starts from 1
1206 // neighbor = neighborlist[i * 4 + j] - 1;
1207 // if (neighbor < 0) {
1208 // face_type[4*i+j] = 5;
1209 // if (cell_type != exterior_cell_type)
1210 // printf("there might be an error with tet attibutes ....\n");
1211 // }
1212 // else {
1213 // neightype = (int)tetrahedronattributelist[neighbor * numberoftetrahedronattributes];
1214 //
1215 // if (ActiveSite != NULL){
1216 // if (neightype != cell_type) {
1217 // active = 1;
1218 // for (k = 0; k < numberofcorners; k++) {
1219 // if (k != j && ActiveSite[tetra_node[k]] == 0)
1220 // active = 0;
1221 // }
1222 // if (active)
1223 // face_type[4*i+j] = ActiveSite[tetra_node[(j+1)%numberofcorners]];
1224 // else
1225 // face_type[4*i+j] = 4;
1226 // }
1227 // }
1228 // }
1229 // }
1230 //
1231 // if (exterior_cell_type == 1) {
1232 // if (cell_type >= 2) cell_type = 1;
1233 // else cell_type = 2;
1234 // }
1235 //
1236 // else if (exterior_cell_type == 2) {
1237 // if (cell_type >= 3) cell_type = 2;
1238 // }
1239 // else {
1240 // printf("there are extra materials inside molecular surfaces\n");
1241 // exit(0);
1242 // }
1243 //
1244 // fprintf(fout,"%10d 0 %5d %5d %5d %5d %5d %10d %10d %10d %10d \n", i, 0,
1245 // face_type[4*i], face_type[4*i+1], face_type[4*i+2], face_type[4*i+3],
1246 // tetra_node[0], tetra_node[1], tetra_node[2], tetra_node[3]);
1247 // }
1248 //
1249 // fprintf(fout,"]; \n");
1250 // fprintf(fout,"mcsf_end=1; \n");
1251 // fclose(fout);
1252 //
1253 // // Release memory
1254 // delete[] node_attr;
1255 // delete[] face_type;
1256 //
1257 //}
1258
1259 /*
1260 * ***************************************************************************
1261 * Routine: write_pdb_gem
1262 *
1263 * Author: Zeyun Yu (zeyun.yu@gmail.com)
1264 *
1265 * Purpose: Write tetrahedral meshes into "Gem_mesh", which can be readily
1266 * converted to FETK-supported format
1267 * ***************************************************************************
1268 */
1269
GemMesh_fromPdb(tetgenio * out,float radius,float centerx,float centery,float centerz,char * ActiveSite,int output_flag)1270 GemMesh* GemMesh_fromPdb(tetgenio* out, float radius, float centerx,
1271 float centery, float centerz, char *ActiveSite,
1272 int output_flag)
1273 {
1274 GemMesh* Gem_mesh;
1275 long i, j, k;
1276 int *NodeAttr;
1277 long NodeIndex;
1278 int *FaceType, TetraNode[4];
1279 float x,y,z;
1280 int type;
1281 float dist;
1282 long a,b,c,d;
1283 float ax,ay,az;
1284 float bx,by,bz;
1285 float cx,cy,cz;
1286 float dx,dy,dz;
1287 int active, ExType;
1288 long neighbor, neightype;
1289 long *map_w_o,*map_w_i;
1290 long out_pnt_num,out_tet_num;
1291 long in_pnt_num,in_tet_num;
1292
1293
1294 // initialization and memory allocation
1295 Gem_mesh = (GemMesh *)malloc(sizeof(GemMesh));
1296 Gem_mesh->dim = 3;
1297 Gem_mesh->dimii = 3;
1298 Gem_mesh->vv = (FETK_VX *)malloc(sizeof(FETK_VX)*out->numberofpoints);
1299 Gem_mesh->ss = (FETK_SS *)malloc(sizeof(FETK_SS)*out->numberoftetrahedra);
1300
1301
1302 // Check the type of exterior tetrahedra
1303 // In Tetgen, it could be radmonly 1 or 2 (or even other integers if more materials are present)
1304 ExType = -1;
1305 for (i = 0; i < out->numberoftetrahedra; i++) {
1306 type = (int)out->tetrahedronattributelist[i * out->numberoftetrahedronattributes];
1307 for (j = 0; j < out->numberofcorners; j++) {
1308 neighbor = out->neighborlist[i * 4 + j] - 1;
1309 if (neighbor < 0) {
1310 ExType = type;
1311 break;
1312 }
1313 }
1314 if (ExType >= 0)
1315 break;
1316 }
1317 printf("ExType = %d \n",ExType);
1318
1319 // assign node attributes
1320 NodeAttr = new int[out->numberofpoints];
1321 for (i = 0; i < out->numberofpoints; i++) {
1322 if (out->pointmarkerlist[i]) {
1323 x = out->pointlist[i * 3];
1324 y = out->pointlist[i * 3 + 1];
1325 z = out->pointlist[i * 3 + 2];
1326 dist = sqrt((x-centerx)*(x-centerx)+(y-centery)*(y-centery)+(z-centerz)*(z-centerz));
1327 if (dist > 0.99*radius)
1328 NodeAttr[i] = 4;
1329 else
1330 NodeAttr[i] = 3;
1331 }
1332 else
1333 NodeAttr[i] = 0;
1334 }
1335
1336 for (i = 0; i < out->numberoftetrahedra; i++) {
1337 type = (int)out->tetrahedronattributelist[i * out->numberoftetrahedronattributes];
1338 for (j = 0; j < out->numberofcorners; j++) {
1339 // Note that index in tetrahedronlist starts from 1 !!!!!!!!
1340 NodeIndex = out->tetrahedronlist[i * out->numberofcorners + j] - 1;
1341
1342 if (NodeAttr[NodeIndex] == 0) {
1343 if (ExType == 1) {
1344 if (type >= 3) type = 2;
1345 NodeAttr[NodeIndex] = ((type==1) ? (2):(1));
1346 }
1347 else if (ExType == 2) {
1348 if (type >= 3) type = 1;
1349 NodeAttr[NodeIndex] = type;
1350 }
1351 else {
1352 printf("there are extra materials inside molecular surfaces\n");
1353 exit(0);
1354 }
1355 }
1356 }
1357 // check tetra orientation
1358 a = out->tetrahedronlist[i * out->numberofcorners]-1;
1359 b = out->tetrahedronlist[i * out->numberofcorners + 1]-1;
1360 c = out->tetrahedronlist[i * out->numberofcorners + 2]-1;
1361 d = out->tetrahedronlist[i * out->numberofcorners + 3]-1;
1362 ax = out->pointlist[a * 3];
1363 ay = out->pointlist[a * 3 + 1];
1364 az = out->pointlist[a * 3 + 2];
1365 bx = out->pointlist[b * 3];
1366 by = out->pointlist[b * 3 + 1];
1367 bz = out->pointlist[b * 3 + 2];
1368 cx = out->pointlist[c * 3];
1369 cy = out->pointlist[c * 3 + 1];
1370 cz = out->pointlist[c * 3 + 2];
1371 dx = out->pointlist[d * 3];
1372 dy = out->pointlist[d * 3 + 1];
1373 dz = out->pointlist[d * 3 + 2];
1374 bx -= ax;
1375 by -= ay;
1376 bz -= az;
1377 cx -= ax;
1378 cy -= ay;
1379 cz -= az;
1380 dx -= ax;
1381 dy -= ay;
1382 dz -= az;
1383 ax = by*cz-bz*cy;
1384 ay = bz*cx-bx*cz;
1385 az = bx*cy-by*cx;
1386 if (ax*dx+ay*dy+az*dz < 0) {
1387 out->tetrahedronlist[i * out->numberofcorners + 1] = c;
1388 out->tetrahedronlist[i * out->numberofcorners + 2] = b;
1389 k = out->neighborlist[i * 4 + 1];
1390 out->neighborlist[i * 4 + 1] = out->neighborlist[i * 4 + 2];
1391 out->neighborlist[i * 4 + 2] = k;
1392 }
1393 }
1394
1395
1396 // output the whole tetra mesh
1397 if (output_flag == 1) {
1398 Gem_mesh->vertices = out->numberofpoints;
1399 Gem_mesh->simplices = out->numberoftetrahedra;
1400 }
1401 out_pnt_num = 0;
1402 in_pnt_num = 0;
1403 for (i = 0; i < out->numberofpoints; i++) {
1404 if (output_flag == 1) {
1405 Gem_mesh->vv[i].id = i;
1406 Gem_mesh->vv[i].chrt = (NodeAttr[i]==4)?(2):(NodeAttr[i]);
1407 Gem_mesh->vv[i].x = out->pointlist[i * 3];
1408 Gem_mesh->vv[i].y = out->pointlist[i * 3 + 1];
1409 Gem_mesh->vv[i].z = out->pointlist[i * 3 + 2];
1410 }
1411 if (NodeAttr[i] != 1)
1412 out_pnt_num++;
1413 if (NodeAttr[i] == 1 || NodeAttr[i] == 3)
1414 in_pnt_num++;
1415 }
1416
1417 out_tet_num = 0;
1418 in_tet_num = 0;
1419 FaceType = new int[out->numberoftetrahedra*4];
1420 for (i = 0; i < out->numberoftetrahedra; i++) {
1421 for (j = 0; j < out->numberofcorners; j++)
1422 // Note that index in tetrahedronlist starts from 1 !!!!!!!!
1423 TetraNode[j] = out->tetrahedronlist[i * out->numberofcorners + j] - 1;
1424
1425 type = (int)out->tetrahedronattributelist[i * out->numberoftetrahedronattributes];
1426 if (type == ExType)
1427 out_tet_num++;
1428 else
1429 in_tet_num++;
1430
1431 // Check face type
1432 for (j = 0; j < out->numberofcorners; j++) {
1433 FaceType[4*i+j] = 0;
1434
1435 // Note that index in neighborlist starts from 1 !!!!!!!!
1436 neighbor = out->neighborlist[i * 4 + j] - 1;
1437 if (neighbor < 0) {
1438 FaceType[4*i+j] = 5;
1439 if (type != ExType)
1440 printf("there might be an error with tet attibutes ....\n");
1441 }
1442 else {
1443 neightype = (int)out->tetrahedronattributelist[neighbor * out->numberoftetrahedronattributes];
1444 if (ActiveSite != NULL){
1445 if (neightype != type) {
1446 active = 1;
1447 for (k = 0; k < out->numberofcorners; k++) {
1448 if (k != j && ActiveSite[TetraNode[k]] == 0)
1449 active = 0;
1450 }
1451 if (active)
1452 FaceType[4*i+j] = ActiveSite[TetraNode[(j+1)%out->numberofcorners]];
1453 else
1454 FaceType[4*i+j] = 4;
1455 }
1456 }
1457 }
1458 }
1459 if (ExType == 1) {
1460 if (type >= 2) type = 1;
1461 else type = 2;
1462 }
1463 else if (ExType == 2) {
1464 if (type >= 3) type = 2;
1465 }
1466 else {
1467 printf("there are extra materials inside molecular surfaces\n");
1468 exit(0);
1469 }
1470 if (output_flag == 1) {
1471 Gem_mesh->ss[i].id = i;
1472 Gem_mesh->ss[i].grp = 0;
1473 Gem_mesh->ss[i].mat = type;
1474 Gem_mesh->ss[i].fa = (FaceType[4*i]==5)?(5):(0);
1475 Gem_mesh->ss[i].fb = (FaceType[4*i+1]==5)?(5):(0);
1476 Gem_mesh->ss[i].fc = (FaceType[4*i+2]==5)?(5):(0);
1477 Gem_mesh->ss[i].fd = (FaceType[4*i+3]==5)?(5):(0);
1478 Gem_mesh->ss[i].na = TetraNode[0];
1479 Gem_mesh->ss[i].nb = TetraNode[1];
1480 Gem_mesh->ss[i].nc = TetraNode[2];
1481 Gem_mesh->ss[i].nd = TetraNode[3];
1482 }
1483 }
1484
1485
1486 // output the inner tetra mesh
1487 if (output_flag == 2) {
1488 Gem_mesh->vertices = in_pnt_num;
1489 Gem_mesh->simplices = in_tet_num;
1490 }
1491 j = 0;
1492 map_w_i = new long[out->numberofpoints];
1493 for (i = 0; i < out->numberofpoints; i++) {
1494 map_w_i[i] = j;
1495 if (NodeAttr[i] == 1 || NodeAttr[i] == 3) {
1496 if (output_flag == 2) {
1497 Gem_mesh->vv[j].id = j;
1498 Gem_mesh->vv[j].chrt = 0;
1499 Gem_mesh->vv[j].x = out->pointlist[i * 3];
1500 Gem_mesh->vv[j].y = out->pointlist[i * 3 + 1];
1501 Gem_mesh->vv[j].z = out->pointlist[i * 3 + 2];
1502 }
1503 j++;
1504 }
1505 }
1506 k = 0;
1507 for (i = 0; i < out->numberoftetrahedra; i++) {
1508 type = (int)out->tetrahedronattributelist[i * out->numberoftetrahedronattributes];
1509 if (type != ExType) {
1510 for (j = 0; j < out->numberofcorners; j++)
1511 TetraNode[j] = map_w_i[out->tetrahedronlist[i * out->numberofcorners + j] - 1];
1512 if (output_flag == 2) {
1513 Gem_mesh->ss[k].id = k;
1514 Gem_mesh->ss[k].grp = 0;
1515 Gem_mesh->ss[k].mat = 0;
1516 Gem_mesh->ss[k].fa = FaceType[i*4];
1517 Gem_mesh->ss[k].fb = FaceType[i*4+1];
1518 Gem_mesh->ss[k].fc = FaceType[i*4+2];
1519 Gem_mesh->ss[k].fd = FaceType[i*4+3];
1520 Gem_mesh->ss[k].na = TetraNode[0];
1521 Gem_mesh->ss[k].nb = TetraNode[1];
1522 Gem_mesh->ss[k].nc = TetraNode[2];
1523 Gem_mesh->ss[k].nd = TetraNode[3];
1524 }
1525 k++;
1526 }
1527 }
1528
1529
1530 // output the outer tetra mesh
1531 if (output_flag == 3) {
1532 Gem_mesh->vertices = out_pnt_num;
1533 Gem_mesh->simplices = out_tet_num;
1534 }
1535 j = 0;
1536 map_w_o = new long[out->numberofpoints];
1537 for (i = 0; i < out->numberofpoints; i++) {
1538 map_w_o[i] = j;
1539 if (NodeAttr[i] != 1) {
1540 if (output_flag == 3) {
1541 Gem_mesh->vv[j].id = j;
1542 Gem_mesh->vv[j].chrt = 0;
1543 Gem_mesh->vv[j].x = out->pointlist[i * 3];
1544 Gem_mesh->vv[j].y = out->pointlist[i * 3 + 1];
1545 Gem_mesh->vv[j].z = out->pointlist[i * 3 + 2];
1546 }
1547 j++;
1548 }
1549 }
1550 k = 0;
1551 for (i = 0; i < out->numberoftetrahedra; i++) {
1552 type = (int)out->tetrahedronattributelist[i * out->numberoftetrahedronattributes];
1553 if (type == ExType) {
1554 for (j = 0; j < out->numberofcorners; j++)
1555 TetraNode[j] = map_w_o[out->tetrahedronlist[i * out->numberofcorners + j] - 1];
1556 if (output_flag == 3) {
1557 Gem_mesh->ss[k].id = k;
1558 Gem_mesh->ss[k].grp = 0;
1559 Gem_mesh->ss[k].mat = 0;
1560 Gem_mesh->ss[k].fa = FaceType[i*4];
1561 Gem_mesh->ss[k].fb = FaceType[i*4+1];
1562 Gem_mesh->ss[k].fc = FaceType[i*4+2];
1563 Gem_mesh->ss[k].fd = FaceType[i*4+3];
1564 Gem_mesh->ss[k].na = TetraNode[0];
1565 Gem_mesh->ss[k].nb = TetraNode[1];
1566 Gem_mesh->ss[k].nc = TetraNode[2];
1567 Gem_mesh->ss[k].nd = TetraNode[3];
1568 }
1569 k++;
1570 }
1571 }
1572
1573 free(NodeAttr);
1574 free(FaceType);
1575 free(map_w_i);
1576 free(map_w_o);
1577
1578 return Gem_mesh;
1579 }
1580
1581
1582 /*
1583 * ***************************************************************************
1584 * Routine: write_pdb_mcsf
1585 *
1586 * Author: Zeyun Yu (zeyun.yu@gmail.com)
1587 *
1588 * Purpose: Write surface/tetrahedral meshes into disks (MCSF/OFF formats)
1589 * ***************************************************************************
1590 */
1591 // FIXME: Depricated. Use write to GemMesh and then to file
1592 //void tetgenio::write_pdb_mcsf(char *inputfilename, float radius, float centerx,
1593 // float centery, float centerz, char *ActiveSite,
1594 // int output_flag)
1595 //{
1596 // long i, j, k;
1597 // FILE *fout;
1598 // int *NodeAttr;
1599 // long NodeIndex;
1600 // int *FaceType, TetraNode[4];
1601 // float x,y,z;
1602 // char file_name[256];
1603 // int type;
1604 // float dist;
1605 // long a,b,c,d;
1606 // float ax,ay,az;
1607 // float bx,by,bz;
1608 // float cx,cy,cz;
1609 // float dx,dy,dz;
1610 // int active, ExType;
1611 // long neighbor, neightype;
1612 // long out_pnt_num,out_tet_num;
1613 // long in_pnt_num,in_tet_num;
1614 // long bem_pnt_num,bem_tri_num;
1615 // long *map_w_o,*map_w_i,*map_w_t;
1616 //
1617 //
1618 //
1619 // // Check the type of exterior tetrahedra
1620 // // In Tetgen, it could be radmonly 1 or 2 (or even other integers if more materials are present)
1621 // ExType = -1;
1622 // for (i = 0; i < numberoftetrahedra; i++) {
1623 // type = (int)tetrahedronattributelist[i * numberoftetrahedronattributes];
1624 // for (j = 0; j < numberofcorners; j++) {
1625 // neighbor = neighborlist[i * 4 + j] - 1;
1626 // if (neighbor < 0) {
1627 // ExType = type;
1628 // break;
1629 // }
1630 // }
1631 // if (ExType >= 0)
1632 // break;
1633 // }
1634 // printf("ExType = %d \n",ExType);
1635 //
1636 // // assign node attributes
1637 // bem_pnt_num = 0;
1638 // NodeAttr = new int[numberofpoints];
1639 // for (i = 0; i < numberofpoints; i++) {
1640 // if (pointmarkerlist[i]) {
1641 // x = pointlist[i * 3];
1642 // y = pointlist[i * 3 + 1];
1643 // z = pointlist[i * 3 + 2];
1644 // dist = sqrt((x-centerx)*(x-centerx)+(y-centery)*(y-centery)+(z-centerz)*(z-centerz));
1645 // if (dist > 0.99*radius)
1646 // NodeAttr[i] = 4;
1647 // else {
1648 // NodeAttr[i] = 3;
1649 // bem_pnt_num++;
1650 // }
1651 // }
1652 // else
1653 // NodeAttr[i] = 0;
1654 // }
1655 //
1656 // for (i = 0; i < numberoftetrahedra; i++) {
1657 // type = (int)tetrahedronattributelist[i * numberoftetrahedronattributes];
1658 // for (j = 0; j < numberofcorners; j++) {
1659 // // Note that index in tetrahedronlist starts from 1 !!!!!!!!
1660 // NodeIndex = tetrahedronlist[i * numberofcorners + j] - 1;
1661 //
1662 // if (NodeAttr[NodeIndex] == 0) {
1663 // if (ExType == 1) {
1664 // if (type >= 3) type = 2;
1665 // NodeAttr[NodeIndex] = ((type==1) ? (2):(1));
1666 // }
1667 // else if (ExType == 2) {
1668 // if (type >= 3) type = 1;
1669 // NodeAttr[NodeIndex] = type;
1670 // }
1671 // else {
1672 // printf("there are extra materials inside molecular surfaces\n");
1673 // exit(0);
1674 // }
1675 // }
1676 // }
1677 // // check tetra orientation
1678 // a = tetrahedronlist[i * numberofcorners]-1;
1679 // b = tetrahedronlist[i * numberofcorners + 1]-1;
1680 // c = tetrahedronlist[i * numberofcorners + 2]-1;
1681 // d = tetrahedronlist[i * numberofcorners + 3]-1;
1682 // ax = pointlist[a * 3];
1683 // ay = pointlist[a * 3 + 1];
1684 // az = pointlist[a * 3 + 2];
1685 // bx = pointlist[b * 3];
1686 // by = pointlist[b * 3 + 1];
1687 // bz = pointlist[b * 3 + 2];
1688 // cx = pointlist[c * 3];
1689 // cy = pointlist[c * 3 + 1];
1690 // cz = pointlist[c * 3 + 2];
1691 // dx = pointlist[d * 3];
1692 // dy = pointlist[d * 3 + 1];
1693 // dz = pointlist[d * 3 + 2];
1694 // bx -= ax;
1695 // by -= ay;
1696 // bz -= az;
1697 // cx -= ax;
1698 // cy -= ay;
1699 // cz -= az;
1700 // dx -= ax;
1701 // dy -= ay;
1702 // dz -= az;
1703 // ax = by*cz-bz*cy;
1704 // ay = bz*cx-bx*cz;
1705 // az = bx*cy-by*cx;
1706 // if (ax*dx+ay*dy+az*dz < 0) {
1707 // tetrahedronlist[i * numberofcorners + 1] = c;
1708 // tetrahedronlist[i * numberofcorners + 2] = b;
1709 // k = neighborlist[i * 4 + 1];
1710 // neighborlist[i * 4 + 1] = neighborlist[i * 4 + 2];
1711 // neighborlist[i * 4 + 2] = k;
1712 // }
1713 // }
1714 //
1715 //
1716 // // output the whole tetra mesh
1717 // if (output_flag == 1) {
1718 // sprintf(file_name, "%s.output.all.m", inputfilename);
1719 // if ((fout=fopen(file_name, "wb"))==NULL){
1720 // printf("write error...\n");
1721 // exit(0);
1722 // };
1723 // fprintf(fout,"mcsf_begin=1; \n \n");
1724 // fprintf(fout," dim = 3; \n");
1725 // fprintf(fout," dimii = 3; \n");
1726 // fprintf(fout," vertices = %d;\n", numberofpoints);
1727 // fprintf(fout,"simplices = %d;\n", numberoftetrahedra);
1728 // fprintf(fout,"vert=[ \n");
1729 // fprintf(fout,"%%------------------------------------------------------------------------------ \n");
1730 // fprintf(fout,"%% Node-ID Chrt X-Coordinate Y-coordinate Z-coordinate \n");
1731 // fprintf(fout,"%%--------- ---- ---------------- ---------------- ---------------- \n");
1732 // }
1733 //
1734 // out_pnt_num = 0;
1735 // in_pnt_num = 0;
1736 // for (i = 0; i < numberofpoints; i++) {
1737 // if (output_flag == 1)
1738 // fprintf(fout, "%10ld %d %17.10e %17.10e %17.10e \n",i,((NodeAttr[i]==4)?(2):(NodeAttr[i])),
1739 // pointlist[i * 3],pointlist[i * 3 + 1],pointlist[i * 3 + 2]);
1740 // if (NodeAttr[i] != 1)
1741 // out_pnt_num++;
1742 // if (NodeAttr[i] == 1 || NodeAttr[i] == 3)
1743 // in_pnt_num++;
1744 // }
1745 //
1746 // if (output_flag == 1) {
1747 // fprintf(fout,"]; \n");
1748 // fprintf(fout,"simp=[ \n");
1749 // fprintf(fout,"%%------------------------------------------------------------------------------------------ \n");
1750 // fprintf(fout,"%% Simp-ID Grp Mat Face-Types Vertex-Numbers \n");
1751 // fprintf(fout,"%%--------- --- --- --------------------- ------------------------------------------- \n");
1752 // }
1753 //
1754 // out_tet_num = 0;
1755 // in_tet_num = 0;
1756 // bem_tri_num = 0;
1757 // FaceType = new int[numberoftetrahedra*4];
1758 // for (i = 0; i < numberoftetrahedra; i++) {
1759 // for (j = 0; j < numberofcorners; j++)
1760 // // Note that index in tetrahedronlist starts from 1 !!!!!!!!
1761 // TetraNode[j] = tetrahedronlist[i * numberofcorners + j] - 1;
1762 //
1763 // type = (int)tetrahedronattributelist[i * numberoftetrahedronattributes];
1764 // if (type == ExType)
1765 // out_tet_num++;
1766 // else
1767 // in_tet_num++;
1768 //
1769 // // Check face type
1770 // for (j = 0; j < numberofcorners; j++) {
1771 // FaceType[4*i+j] = 0;
1772 //
1773 // // Note that index in neighborlist starts from 1 !!!!!!!!
1774 // neighbor = neighborlist[i * 4 + j] - 1;
1775 // if (neighbor < 0) {
1776 // FaceType[4*i+j] = 5;
1777 // if (type != ExType)
1778 // printf("there might be an error with tet attibutes ....\n");
1779 // }
1780 // else {
1781 // neightype = (int)tetrahedronattributelist[neighbor * numberoftetrahedronattributes];
1782 //
1783 // if (ActiveSite != NULL){
1784 // if (neightype != type) {
1785 // active = 1;
1786 // for (k = 0; k < numberofcorners; k++) {
1787 // if (k != j && ActiveSite[TetraNode[k]] == 0)
1788 // active = 0;
1789 // }
1790 // if (active)
1791 // FaceType[4*i+j] = ActiveSite[TetraNode[(j+1)%numberofcorners]];
1792 // else
1793 // FaceType[4*i+j] = 4;
1794 // }
1795 // }
1796 // }
1797 // if (type == ExType && (FaceType[4*i+j] != 0 && FaceType[4*i+j] != 5))
1798 // bem_tri_num++;
1799 //
1800 // }
1801 // if (ExType == 1) {
1802 // if (type >= 2) type = 1;
1803 // else type = 2;
1804 // }
1805 // else if (ExType == 2) {
1806 // if (type >= 3) type = 2;
1807 // }
1808 // else {
1809 // printf("there are extra materials inside molecular surfaces\n");
1810 // exit(0);
1811 // }
1812 // if (output_flag == 1)
1813 // fprintf(fout,"%10ld 0 %5d %5d %5d %5d %5d %10d %10d %10d %10d \n",i,type,((FaceType[4*i]==5)?(5):(0)),
1814 // ((FaceType[4*i+1]==5)?(5):(0)),((FaceType[4*i+2]==5)?(5):(0)),((FaceType[4*i+3]==5)?(5):(0)),
1815 // TetraNode[0],TetraNode[1],TetraNode[2],TetraNode[3]);
1816 // }
1817 // if (output_flag == 1) {
1818 // fprintf(fout,"]; \n");
1819 // fprintf(fout,"mcsf_end=1; \n");
1820 // fclose(fout);
1821 // }
1822 //
1823 //
1824 // // output the inner tetra mesh
1825 // if (output_flag == 2) {
1826 // sprintf(file_name, "%s.output.in.m", inputfilename);
1827 // if ((fout=fopen(file_name, "wb"))==NULL){
1828 // printf("write error...\n");
1829 // exit(0);
1830 // };
1831 // fprintf(fout,"mcsf_begin=1; \n \n");
1832 // fprintf(fout," dim = 3; \n");
1833 // fprintf(fout," dimii = 3; \n");
1834 // fprintf(fout," vertices = %ld;\n", in_pnt_num);
1835 // fprintf(fout,"simplices = %ld;\n", in_tet_num);
1836 // fprintf(fout,"vert=[ \n");
1837 // fprintf(fout,"%%------------------------------------------------------------------------------------ \n");
1838 // fprintf(fout,"%% Node-ID Chrt X-Coordinate Y-coordinate Z-coordinate \n");
1839 // fprintf(fout,"%%--------- ---------- ---------------- ---------------- ---------------- \n");
1840 // }
1841 // j = 0;
1842 // map_w_i = new long[numberofpoints];
1843 // for (i = 0; i < numberofpoints; i++) {
1844 // map_w_i[i] = j;
1845 // if (NodeAttr[i] == 1 || NodeAttr[i] == 3) {
1846 // if (output_flag == 2)
1847 // fprintf(fout, "%10ld %10d %17.10e %17.10e %17.10e \n",j,0,
1848 // pointlist[i * 3],pointlist[i * 3 + 1],pointlist[i * 3 + 2]);
1849 // j++;
1850 // }
1851 // }
1852 // if (output_flag == 2) {
1853 // fprintf(fout,"]; \n");
1854 // fprintf(fout,"simp=[ \n");
1855 // fprintf(fout,"%%--------------------------------------------------------------------------------------- \n");
1856 // fprintf(fout,"%% Simp-ID Grp Mat Face-Types Vertex-Numbers \n");
1857 // fprintf(fout,"%%--------- --- --- --------------------- ------------------------------------------- \n");
1858 // }
1859 // k = 0;
1860 // for (i = 0; i < numberoftetrahedra; i++) {
1861 // type = (int)tetrahedronattributelist[i * numberoftetrahedronattributes];
1862 // if (type != ExType) {
1863 // for (j = 0; j < numberofcorners; j++)
1864 // TetraNode[j] = map_w_i[tetrahedronlist[i * numberofcorners + j] - 1];
1865 // if (output_flag == 2)
1866 // fprintf(fout,"%10ld 0 0 %5d %5d %5d %5d %10d %10d %10d %10d \n",k,FaceType[i*4],FaceType[i*4+1],
1867 // FaceType[i*4+2],FaceType[i*4+3],TetraNode[0],TetraNode[1],TetraNode[2],TetraNode[3]);
1868 // k++;
1869 // }
1870 // }
1871 // if (output_flag == 2) {
1872 // fprintf(fout,"]; \n");
1873 // fprintf(fout,"mcsf_end=1; \n");
1874 // fclose(fout);
1875 // }
1876 //
1877 //
1878 // // output the outer tetra mesh
1879 // if (output_flag == 3) {
1880 // sprintf(file_name, "%s.output.out.m", inputfilename);
1881 // if ((fout=fopen(file_name, "wb"))==NULL){
1882 // printf("write error...\n");
1883 // exit(0);
1884 // };
1885 // fprintf(fout,"mcsf_begin=1; \n \n");
1886 // fprintf(fout," dim = 3; \n");
1887 // fprintf(fout," dimii = 3; \n");
1888 // fprintf(fout," vertices = %ld;\n", out_pnt_num);
1889 // fprintf(fout,"simplices = %ld;\n", out_tet_num);
1890 // fprintf(fout,"vert=[ \n");
1891 // fprintf(fout,"%%------------------------------------------------------------------------------------ \n");
1892 // fprintf(fout,"%% Node-ID Chrt X-Coordinate Y-coordinate Z-coordinate \n");
1893 // fprintf(fout,"%%--------- ---------- ---------------- ---------------- ---------------- \n");
1894 // }
1895 // j = 0;
1896 // map_w_o = new long[numberofpoints];
1897 // for (i = 0; i < numberofpoints; i++) {
1898 // map_w_o[i] = j;
1899 // if (NodeAttr[i] != 1) {
1900 // if (output_flag == 3)
1901 // fprintf(fout, "%10ld %10d %17.10e %17.10e %17.10e \n",j,0,
1902 // pointlist[i * 3],pointlist[i * 3 + 1],pointlist[i * 3 + 2]);
1903 // j++;
1904 // }
1905 // }
1906 // if (output_flag == 3) {
1907 // fprintf(fout,"]; \n");
1908 // fprintf(fout,"simp=[ \n");
1909 // fprintf(fout,"%%--------------------------------------------------------------------------------------- \n");
1910 // fprintf(fout,"%% Simp-ID Grp Mat Face-Types Vertex-Numbers \n");
1911 // fprintf(fout,"%%--------- --- --- --------------------- ------------------------------------------- \n");
1912 // }
1913 // k = 0;
1914 // for (i = 0; i < numberoftetrahedra; i++) {
1915 // type = (int)tetrahedronattributelist[i * numberoftetrahedronattributes];
1916 // if (type == ExType) {
1917 // for (j = 0; j < numberofcorners; j++)
1918 // TetraNode[j] = map_w_o[tetrahedronlist[i * numberofcorners + j] - 1];
1919 // if (output_flag == 3)
1920 // fprintf(fout,"%10ld 0 0 %5d %5d %5d %5d %10d %10d %10d %10d \n",k,FaceType[i*4],FaceType[i*4+1],
1921 // FaceType[i*4+2],FaceType[i*4+3],TetraNode[0],TetraNode[1],TetraNode[2],TetraNode[3]);
1922 // k++;
1923 // }
1924 // }
1925 // if (output_flag == 3) {
1926 // fprintf(fout,"]; \n");
1927 // fprintf(fout,"mcsf_end=1; \n");
1928 // fclose(fout);
1929 // }
1930
1931
1932
1933 // the following may be useful for boundary element analysis
1934 /*
1935 // output the mapping file from outer domain to whole domain
1936 sprintf(file_name, "%s.output.map.m", inputfilename);
1937 if ((fout=fopen(file_name, "wb"))==NULL){
1938 printf("write error...\n");
1939 exit(0);
1940 };
1941 j = 0;
1942 for (i = 0; i < numberofpoints; i++) {
1943 if (NodeAttr[i] != 1) {
1944 fprintf(fout,"%10d %10d \n",j+1,i);
1945 j++;
1946 }
1947 }
1948 fclose(fout);
1949
1950
1951 // output the molecular surface mesh in .m
1952 sprintf(file_name, "%s.output.bem.m", inputfilename);
1953 if ((fout=fopen(file_name, "wb"))==NULL){
1954 printf("write error...\n");
1955 exit(0);
1956 };
1957 fprintf(fout,"%10d %10d 0 \n",bem_pnt_num,bem_tri_num);
1958 map_w_t = new long[numberofpoints];
1959 j = 0;
1960 for (i = 0; i < numberofpoints; i++) {
1961 map_w_t[i] = j;
1962 if (NodeAttr[i] == 3) {
1963 fprintf(fout,"%13.4f %13.4f %13.4f %10d \n",pointlist[i * 3],pointlist[i * 3 + 1],
1964 pointlist[i * 3 + 2],map_w_o[i]+1);
1965 j++;
1966 }
1967 }
1968 for (i = 0; i < numberoftetrahedra; i++) {
1969 type = (int)tetrahedronattributelist[i * numberoftetrahedronattributes];
1970 if (type == ExType) {
1971 for (j = 0; j < numberofcorners; j++) {
1972 if (FaceType[4*i+j] != 0 && FaceType[4*i+j] != 5) {
1973 a = 0;
1974 for (k = 0; k < numberofcorners; k++) {
1975 if (k != j) {
1976 TetraNode[a] = map_w_t[tetrahedronlist[i * numberofcorners + k] - 1];
1977 a++;
1978 }
1979 }
1980 if (j == 0 || j == 2)
1981 fprintf(fout,"%10d %10d %10d \n",TetraNode[0]+1,TetraNode[1]+1,TetraNode[2]+1);
1982 else if (j == 1 || j == 3)
1983 fprintf(fout,"%10d %10d %10d \n",TetraNode[0]+1,TetraNode[2]+1,TetraNode[1]+1);
1984 }
1985 }
1986 }
1987 }
1988 fclose(fout);
1989
1990
1991 // randomly assign colors to different types of active sites
1992 float red[256],green[256],blue[256];
1993 for (i = 0; i < 256; i++) {
1994 red[i] = (float) rand()/(float)RAND_MAX;
1995 green[i] = (float) rand()/(float)RAND_MAX;
1996 blue[i] = (float) rand()/(float)RAND_MAX;
1997 }
1998 // non-active site color
1999 red[4] = 0.5;
2000 green[4] = 1;
2001 blue[4] = 1;
2002 // active site colors
2003 red[1] = 1;
2004 green[1] = 1;
2005 blue[1] = 0.1;
2006 red[3] = 0.1;
2007 green[3] = 0.1;
2008 blue[3] = 1;
2009 red[7] = 0.1;
2010 green[7] = 1;
2011 blue[7] = 0.1;
2012 red[9] = 1;
2013 green[9] = 0.1;
2014 blue[9] = 0.1;
2015 red[11] = 1;
2016 green[11] = 0.1;
2017 blue[11] = 1;
2018 red[13] = 1;
2019 green[13] = 0.5;
2020 blue[13] = 1;
2021 // output the molecular surface mesh in .off
2022 sprintf(file_name, "%s.output.bem.off", inputfilename);
2023 if ((fout=fopen(file_name, "wb"))==NULL){
2024 printf("write error...\n");
2025 exit(0);
2026 };
2027 fprintf(fout,"OFF\n%10d %10d 0 \n",bem_pnt_num,bem_tri_num);
2028 map_w_t = new long[numberofpoints];
2029 j = 0;
2030 for (i = 0; i < numberofpoints; i++) {
2031 map_w_t[i] = j;
2032 if (NodeAttr[i] == 3) {
2033 fprintf(fout,"%13.4f %13.4f %13.4f \n",pointlist[i * 3],pointlist[i * 3 + 1],
2034 pointlist[i * 3 + 2]);
2035 j++;
2036 }
2037 }
2038 for (i = 0; i < numberoftetrahedra; i++) {
2039 type = (int)tetrahedronattributelist[i * numberoftetrahedronattributes];
2040 if (type == ExType) {
2041 for (j = 0; j < numberofcorners; j++) {
2042 if (FaceType[4*i+j] != 0 && FaceType[4*i+j] != 5) {
2043 a = 0;
2044 for (k = 0; k < numberofcorners; k++) {
2045 if (k != j) {
2046 TetraNode[a] = map_w_t[tetrahedronlist[i * numberofcorners + k] - 1];
2047 a++;
2048 }
2049 }
2050 if (j == 0 || j == 2)
2051 fprintf(fout,"3 %10d %10d %10d %f %f %f \n",TetraNode[0],TetraNode[1],TetraNode[2],
2052 red[FaceType[4*i+j]],green[FaceType[4*i+j]],blue[FaceType[4*i+j]]);
2053 else if (j == 1 || j == 3)
2054 fprintf(fout,"3 %10d %10d %10d %f %f %f \n",TetraNode[0],TetraNode[2],TetraNode[1],
2055 red[FaceType[4*i+j]],green[FaceType[4*i+j]],blue[FaceType[4*i+j]]);
2056 }
2057 }
2058 }
2059 }
2060 fclose(fout);
2061 */
2062 //}
2063
2064 /*
2065 * ***************************************************************************
2066 * Routine: ReadActiveSiteFile
2067 *
2068 * Author: Johan Hake (joha.hake@gmail.com), Zeyun Yu (zeyun.yu@gmail.com)
2069 *
2070 * Purpose: Read a file containing information about active sites
2071 * ***************************************************************************
2072 */
ReadActiveSiteFile(char * input_site,unsigned int & num_spheres,ATOM * & sphere_list,unsigned int * & sphere_markers)2073 void ReadActiveSiteFile(char* input_site,
2074 unsigned int& num_spheres,
2075 ATOM*& sphere_list, unsigned int*& sphere_markers)
2076 // ,float*& area_constraint_list)
2077 {
2078 FILE *fout;
2079 char buf[1024];
2080 int i;
2081 unsigned int marker, c;
2082 float x, y, z, radius, area_constraint;
2083 char* ret_gets;
2084 int ret_scan;
2085
2086 printf("Reading active site file\n");
2087 if ((fout=fopen(input_site, "r"))==NULL){
2088 printf("read error...\n");
2089 exit(0);
2090 };
2091 while ((c=fgetc(fout)) == '#')
2092 ret_gets = fgets(buf, 1024, fout);
2093 ungetc(c, fout);
2094 ret_scan = fscanf(fout,"%d\n", &num_spheres);
2095 sphere_list = (ATOM*)malloc(sizeof(ATOM)*num_spheres);
2096 sphere_markers = (unsigned int*)malloc(sizeof(unsigned int)*num_spheres);
2097 //area_constraint_list = (float*)malloc(sizeof(float)*num_spheres);
2098 for (i = 0; i < num_spheres; i++) {
2099 //fscanf(fout, "%f %f %f %f %d\n", &x, &y, &z, &radius, &marker, &area_constraint);
2100 ret_scan = fscanf(fout, "%f %f %f %f %d\n", &x, &y, &z, &radius, &marker);
2101 sphere_list[i].x = x;
2102 sphere_list[i].y = y;
2103 sphere_list[i].z = z;
2104 sphere_list[i].radius = radius;
2105 sphere_markers[i] = marker;
2106 //area_constraint_list[i] = area_constraint;
2107 }
2108 fclose(fout);
2109 }
2110
2111