/**
# Additional output functions to generate surface values */
#include "geometry.h"
#include "fractions.h"
#if dimension == 1
coord mycs (Point point, scalar c) {
coord n = {1.};
return n;
}
#elif dimension == 2
# include "myc2d.h"
#else // dimension == 3
# include "myc.h"
#endif
void output_ply (struct OutputFacets p)
{
#if defined(_OPENMP)
int num_omp = omp_get_max_threads();
omp_set_num_threads(1);
#endif
scalar c = p.c;
face vector s = p.s;
if (!p.fp) p.fp = stdout;
// print header text
fputs ("ply\n", p.fp);
fputs ("format ascii 1.0\n", p.fp);
int nverts = 0;
int nfacets = 0;
foreach()
if (c[] > 1e-6 && c[] < 1. - 1e-6) {
coord n;
if (!s.x.i)
// compute normal from volume fraction
n = mycs (point, c);
else {
// compute normal from face fractions
double nn = 0.;
foreach_dimension() {
n.x = s.x[] - s.x[1];
nn += fabs(n.x);
}
assert (nn > 0.);
foreach_dimension()
n.x /= nn;
}
double alpha = plane_alpha (c[], n);
coord v[12];
int m = facets (n, alpha, v, 1.);
for (int i = 0; i < m; i++) {
nverts ++;
}
if (m > 0) {
nfacets ++;
}
}
fprintf (p.fp, "element vertex %i\n", nverts);
fputs ("property float x\n", p.fp);
fputs ("property float y\n", p.fp);
fputs ("property float z\n", p.fp);
fprintf (p.fp, "element face %i\n", nfacets);
fputs ("property list uchar int vertex_index\n", p.fp);
fputs ("end_header\n", p.fp);
int facet_num[nfacets];
int ifacet = 0;
int ivert = 0;
foreach()
if (c[] > 1e-6 && c[] < 1. - 1e-6) {
coord n;
if (!s.x.i)
// compute normal from volume fraction
n = mycs (point, c);
else {
// compute normal from face fractions
double nn = 0.;
foreach_dimension() {
n.x = s.x[] - s.x[1];
nn += fabs(n.x);
}
assert (nn > 0.);
foreach_dimension()
n.x /= nn;
}
double alpha = plane_alpha (c[], n);
coord v[12];
int m = facets (n, alpha, v, 1.);
for (int i = 0; i < m; i++) {
fprintf (p.fp, "%g %g %g\n",
x + v[i].x*Delta, y + v[i].y*Delta, z + v[i].z*Delta);
}
if (m > 0) {
facet_num[ifacet] = m;
ifacet ++;
}
}
// print face list
for (ifacet = 0; ifacet < nfacets; ifacet++) {
fprintf (p.fp, "%i ", facet_num[ifacet]);
for (int iv = 0; iv < facet_num[ifacet]; iv ++) {
fprintf (p.fp, "%i ", ivert);
ivert ++;
}
fputc ('\n', p.fp);
}
fflush (p.fp);
#if defined(_OPENMP)
omp_set_num_threads(num_omp);
#endif
}
void output_vtu (struct OutputFacets p)
{
#if defined(_OPENMP)
int num_omp = omp_get_max_threads();
omp_set_num_threads(1);
#endif
scalar c = p.c;
face vector s = p.s;
if (!p.fp) p.fp = stdout;
// print header text
fputs ("\n", p.fp);
fputs ("\n", p.fp);
fputs (" \n", p.fp);
int nverts = 0;
int nfacets = 0;
foreach()
if (c[] > 1e-6 && c[] < 1. - 1e-6) {
coord n;
if (!s.x.i)
// compute normal from volume fraction
n = mycs (point, c);
else {
// compute normal from face fractions
double nn = 0.;
foreach_dimension() {
n.x = s.x[] - s.x[1];
nn += fabs(n.x);
}
assert (nn > 0.);
foreach_dimension()
n.x /= nn;
}
double alpha = plane_alpha (c[], n);
coord v[12];
int m = facets (n, alpha, v, 1.);
for (int i = 0; i < m; i++) {
nverts ++;
}
if (m > 0) {
nfacets ++;
}
}
fprintf (p.fp, " \n", nverts, nfacets);
fputs (" \n", p.fp);
fputs (" \n", p.fp);
int offsets[nfacets];
int ifacet = 0;
int offset = 0;
foreach()
if (c[] > 1e-6 && c[] < 1. - 1e-6) {
coord n;
if (!s.x.i)
// compute normal from volume fraction
n = mycs (point, c);
else {
// compute normal from face fractions
double nn = 0.;
foreach_dimension() {
n.x = s.x[] - s.x[1];
nn += fabs(n.x);
}
assert (nn > 0.);
foreach_dimension()
n.x /= nn;
}
double alpha = plane_alpha (c[], n);
coord v[12];
int m = facets (n, alpha, v, 1.);
for (int i = 0; i < m; i++) {
fprintf (p.fp, "%g %g %g ",
x + v[i].x*Delta, y + v[i].y*Delta, z + v[i].z*Delta);
}
if (m > 0) {
offset += m;
offsets[ifacet] = offset;
ifacet ++;
}
}
fputs (" \n", p.fp);
fputs (" \n", p.fp);
fputs (" \n", p.fp);
fputs (" \n", p.fp);
// print vert numbers
for (int ivert = 0; ivert < nverts; ivert++)
fprintf (p.fp, "%i ", ivert);
fputs (" \n", p.fp);
fputs (" \n", p.fp);
// print offsets
for (ifacet = 0; ifacet < nfacets; ifacet++)
fprintf (p.fp, "%i ", offsets[ifacet]);
fputs (" \n", p.fp);
fputs (" \n", p.fp);
// print cell type list
for (ifacet = 0; ifacet < nfacets; ifacet++)
fprintf (p.fp, "7 ");
fputs (" \n", p.fp);
fputs (" \n", p.fp);
fputs (" \n", p.fp);
fputs (" \n", p.fp);
fputs (" \n", p.fp);
fputs (" \n", p.fp);
fputs (" \n", p.fp);
fputs (" \n", p.fp);
fputs ("\n", p.fp);
fflush (p.fp);
#if defined(_OPENMP)
omp_set_num_threads(num_omp);
#endif
}
struct OutputFacets_scalar {
scalar c;
FILE * fp; // optional: default is stdout
scalar * list; // List of scalar fields to include when writing vtu surface to file
vector * vlist; // List of vector fields to include.
face vector s; // optional: default is none
};
void output_vtu_w_fielddata (struct OutputFacets_scalar p)
{
#if defined(_OPENMP)
int num_omp = omp_get_max_threads();
omp_set_num_threads(1);
#endif
scalar c = p.c;
face vector s = p.s;
if (!p.fp) p.fp = stdout;
// print header text
fputs ("\n", p.fp);
fputs ("\n", p.fp);
fputs ("\t\n", p.fp);
int nverts = 0;
int nfacets = 0;
foreach()
if (c[] > 1e-6 && c[] < 1. - 1e-6) {
coord n;
if (!s.x.i)
// compute normal from volume fraction
n = mycs (point, c);
else {
// compute normal from face fractions
double nn = 0.;
foreach_dimension() {
n.x = s.x[] - s.x[1];
nn += fabs(n.x);
}
assert (nn > 0.);
foreach_dimension()
n.x /= nn;
}
double alpha = plane_alpha (c[], n);
coord v[12];
int m = facets (n, alpha, v, 1.);
for (int i = 0; i < m; i++) {
nverts ++;
}
if (m > 0) {
nfacets ++;
}
}
fprintf (p.fp, "\t\t\n", nverts, nfacets);
// Write list of scalar field values to file
fputs ("\t\t\t \n", p.fp);
for (scalar s in p.list) {
fprintf (p.fp,"\t\t\t\t \n", s.name);
foreach(){
if (c[] > 1e-6 && c[] < 1. - 1e-6) {
fprintf (p.fp, "%g\n", val(s));
}
}
fputs ("\t\t\t\t \n", p.fp);
}
// Write list of vector field values to file
for (vector v in p.vlist) {
fprintf (p.fp,"\t\t\t\t \n", v.x.name);
foreach(){
if (c[] > 1e-6 && c[] < 1. - 1e-6) {
#if dimension == 2
fprintf (p.fp, "%g %g 0.\n", val(v.x), val(v.y));
#endif
#if dimension == 3
fprintf (p.fp, "%g %g %g\n", val(v.x), val(v.y), val(v.z));
#endif
}
}
fputs ("\t\t\t\t \n", p.fp);
}
fputs ("\t\t\t \n", p.fp);
// Write points to file
fputs (" \n", p.fp);
fputs (" \n", p.fp);
int offsets[nfacets];
int ifacet = 0;
int offset = 0;
foreach()
if (c[] > 1e-6 && c[] < 1. - 1e-6) {
coord n;
if (!s.x.i)
// compute normal from volume fraction
n = mycs (point, c);
else {
// compute normal from face fractions
double nn = 0.;
foreach_dimension() {
n.x = s.x[] - s.x[1];
nn += fabs(n.x);
}
assert (nn > 0.);
foreach_dimension()
n.x /= nn;
}
double alpha = plane_alpha (c[], n);
coord v[12];
int m = facets (n, alpha, v, 1.);
for (int i = 0; i < m; i++) {
fprintf (p.fp, "%g %g %g ",
x + v[i].x*Delta, y + v[i].y*Delta, z + v[i].z*Delta);
}
if (m > 0) {
offset += m;
offsets[ifacet] = offset;
ifacet ++;
}
}
fputs (" \n", p.fp);
fputs (" \n", p.fp);
fputs (" \n", p.fp);
fputs (" \n", p.fp);
// print vert numbers
for (int ivert = 0; ivert < nverts; ivert++)
fprintf (p.fp, "%i ", ivert);
fputs (" \n", p.fp);
fputs (" \n", p.fp);
// print offsets
for (ifacet = 0; ifacet < nfacets; ifacet++)
fprintf (p.fp, "%i ", offsets[ifacet]);
fputs (" \n", p.fp);
fputs (" \n", p.fp);
// print cell type list
for (ifacet = 0; ifacet < nfacets; ifacet++)
fprintf (p.fp, "7 ");
fputs (" \n", p.fp);
fputs (" \n", p.fp);
fputs (" \n", p.fp);
fputs (" \n", p.fp);
//fputs (" \n", p.fp);
//fputs (" \n", p.fp);
fputs (" \n", p.fp);
fputs (" \n", p.fp);
fputs ("\n", p.fp);
fflush (p.fp);
#if defined(_OPENMP)
omp_set_num_threads(num_omp);
#endif
}