output_vts.h

/*
This function outputs the entire computational domain from a multilayer solver to a vts file which then can be read using Paraview. Is compatible with multigrid and OpenMP (Thanks to Øystein Lande, whose code this function was developed from.) 
*/
    

void output_vts_ascii_all_layers(FILE* fp, scalar* list, int N)
{
    int nthreads_ = omp_get_num_threads();
    omp_set_num_threads(1); 
    
        // MULTIGRID and preamble. Defining metadata

 fputs("<?xml version=\"1.0\"?>\n<VTKFile type=\"StructuredGrid\" version=\"1.0\" byte_order=\"LittleEndian\" header_type=\"UInt64\">\n", fp);

#if dimension == 1
    fprintf(fp, "\t <StructuredGrid WholeExtent=\"%d %d %d %d 0 0\">\n", 0, N, 0, nl);
    fprintf(fp, "\t\t <Piece Extent=\"%d %d %d %d  0 0\">\n", 0, N, 0, nl);
#endif

#if dimension == 2
    fprintf(fp, "\t <StructuredGrid WholeExtent=\"%d %d %d %d %d %d\">\n", 0, N, 0, N, 0, nl);
    fprintf(fp, "\t\t <Piece Extent=\"%d %d %d %d %d %d\">\n", 0, N, 0, N, 0, nl);
#endif

    // Loop over velocity data and store kinematics in cell vector stucture
    fputs("\t\t\t <CellData Scalars=\"scalars\">\n", fp);
    fprintf(fp, "\t\t\t\t <DataArray type=\"Float64\" NumberOfComponents=\"3\" Name=\"Velocity\" format=\"ascii\">\n");
    for(int i = nl-1; i >= 0; i--){
    foreach() {
#if dimension == 1
        fprintf(fp, "%g %g 0.\n", u.x[0,0,i], w[0,0,i]);
#endif
#if dimension == 2
        fprintf(fp, "%g %g %g\n", u.x[0,0,i], u.y[0,0,i], w[0,0,i]);
#endif
    }
    }
    // Dummy text to get correct amount of layers
    foreach() {
        fprintf(fp, "0 0 0.\n");
    }
    fputs("\t\t\t\t </DataArray>\n", fp);


    // loop over all scalars in scalarlist and store values as cell data
    for (scalar s in list) {
        if (strcmp(s.name, "eta") == 0) {
            fprintf(fp, "\t\t\t\t <DataArray type=\"Float64\" Name=\"%s\" format=\"ascii\">\n", s.name);
            for(int i = nl-1; i >= 0; i--){
                foreach() {
                if (h[] > dry) {
                    fprintf(fp, "%g\n", s[0,0,i]);
                }
                else {
                    fprintf(fp, "nan\n");
                }
            }
        }
        fputs("\t\t\t\t </DataArray>\n", fp);
        }
        else {
            fprintf(fp, "\t\t\t\t <DataArray type=\"Float64\" Name=\"%s\" format=\"ascii\">\n", s.name);
            for(int i = nl-1; i >= 0; i--){
            foreach() {
                fprintf(fp, "%g\n", s[0,0,i]);
            }
            }
            foreach(){
                fprintf(fp, "0\n");
            }
            fputs("\t\t\t\t </DataArray>\n", fp);
        }   
    }   

    fputs("\t\t\t </CellData>\n", fp);
    // Coordinates, maps height using layer thickness 
    fputs("\t\t\t <Points>\n", fp);
#if dimension == 1
    fputs("\t\t\t\t <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">\n", fp);
        double* zcorr = (double *)malloc((N+1)*sizeof(double));
        for (int j = 0; j <= N; j++){
            zcorr[j] = 0;
        }
        int j;
        for(int i = nl-1; i >= 0; i--){
            j = 0;
            foreach_vertex(serial) {
                fprintf(fp, "%12.4f %12.4f 0.\n", x,zcorr[j],0 );
                zcorr[j] = zcorr[j] - h[0,0,i];
                if (h[0,0i] < 1e-5){
                  zcorr[j] = h[-1,-1,i];  //Correction for vertices on domain edge where h[0,0,i] \approx 0.
                }
                j++;
            }
        }
        foreach_vertex(serial){
                fprintf(fp, "%12.4f %12.4f 0.\n", x, zcorr);
        }
        free(zcorr);
#endif
#if dimension == 2
    fputs("\t\t\t\t <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">\n", fp);
        
        double* zcorr = (double*)malloc((N+1)*(N+1)*sizeof(double));
        for (int j = 0; j < (N+1)*(N+1); j++){
            zcorr[j] = 0;
        }
        int j, k;
        for(int i = nl-1; i >= 0; i--){
            j = 0;
            foreach_vertex(serial) {
                fprintf(fp, "%12.4f %12.4f %12.4f\n", x,y, zcorr[j]);
                zcorr[j] = zcorr[j] - h[0,0,i];
                if (h[0,0,i] < 1e-3){
                    zcorr[j] = h[-1,-1,i]; //Correction for vertices on domain edge where h[0,0,i] \approx 0.
                }
                j++;
         }
        j = 0;
        foreach_vertex(serial){
                fprintf(fp, "%12.4f %12.4f %12.4f\n", x,y, zcorr[j]);
                j++;
        }
        }
        free(zcorr);
#endif

    fputs("\t\t\t\t </DataArray>\n", fp);
    fputs("\t\t\t </Points>\n", fp);
    fputs("\t\t </Piece>\n", fp);

    // write time value
    fprintf(fp, "\t\t <FieldData> \n");
    fprintf(fp, "\t\t\t <DataArray type = \"Float64\" Name = \"%s\" NumberOfTuples = \"1\" format = \"ascii\" RangeMin = \"%.3f\" RangeMax = \"%.3f\"> \n", "TimeValue", t + dt_start, t + dt_start);
    fprintf(fp, "\t\t\t %.3f \n", t + dt_start);
    fprintf(fp, "\t\t\t </DataArray > \n");
    fprintf(fp, "\t\t </FieldData> \n");
    fputs("\t </StructuredGrid>\n", fp);
    fputs("</VTKFile>\n", fp);
    fflush(fp);

    omp_set_num_threads(nthreads_);

}