fix minor bugs when writing out avs filev0.2

8 files changed, 384 insertions, 12 deletions

diff --git a/Makefile b/Makefile
index 6840954..597d169 100644
--- a/Makefile
+++ b/Makefile
@@ -19,7 +19,7 @@ CTAGS=uctags
 TRIANGLE_DIR=/Users/IOzgen/Documents/work/coding/tools/triangle
 
 CFLAGS=-g -Wall -Wextra -pedantic -std=c99
-LFLAGS=-I$(TRIANGLE_DIR)
+LFLAGS=-I$(TRIANGLE_DIR) -fopenmp
 CPU=
 RUNFLAGS=
 INCLUDES=
diff --git a/header.h b/header.h
index cfb4162..172649c 100644
--- a/header.h
+++ b/header.h
@@ -106,14 +106,19 @@ int matmul(double* matrix1, int nrow1, int ncol1, double* matrix2, int nrow2, in
   }
   
   double sum = 0.0;
-
-  for (int i = 0; i < nrow1; i ++) {
-    for (int j = 0; j < ncol2; j ++) {
-      for (int k = 0; k < ncol1; k ++) {
-	sum += matrix1[k + i * ncol1] * matrix2[j + k * ncol2];
+  int i;
+
+#pragma omp parallel
+  {
+#pragma omp for private(i) reduction(+:sum)
+    for (i = 0; i < nrow1; i ++) {
+      for (int j = 0; j < ncol2; j ++) {
+	for (int k = 0; k < ncol1; k ++) {
+	  sum += matrix1[k + i * ncol1] * matrix2[j + k * ncol2];
+	}
+	result[j + ncol2 * i] = sum;
+	sum = 0.0;
       }
-      result[j + ncol2 * i] = sum;
-      sum = 0.0;
     }
   }
 
@@ -236,7 +241,7 @@ void write_avs(struct triangulateio *io, int counter, double *elevations, int nc
   FILE *fp;
   fp = fopen(fname, "w");
 
-  fprintf(fp, "%d %d %d %d %d\n", io->numberofpoints, io->numberoftriangles, 0, 1, 0);
+  fprintf(fp, "%d %d %d %d %d\n", io->numberofpoints, io->numberoftriangles, 0, 0, 0);
 
   for (int i = 0; i < io->numberofpoints; ++ i) {
 
@@ -277,12 +282,12 @@ void write_avs(struct triangulateio *io, int counter, double *elevations, int nc
       
     }
     
-    fprintf(fp, "%d %f %f %f\n", i+1, io->pointlist[2*i], io->pointlist[2*i+1], z);
+    fprintf(fp, "%08d %f %f %f\n", i+1, io->pointlist[2*i], io->pointlist[2*i+1], z);
 
   }
 
   for (int i = 0; i < io->numberoftriangles; ++ i) {
-    fprintf(fp, "%d %d tri %d %d %d\n", i+1, 1, io->trianglelist[3*i], io->trianglelist[3*i+1], io->trianglelist[3*i+2]);
+    fprintf(fp, "%08d %8d tri %8d %8d %8d\n", i+1, 1, io->trianglelist[3*i], io->trianglelist[3*i+1], io->trianglelist[3*i+2]);
   }
 
   fclose(fp);
diff --git a/main.c b/main.c
index 82ac8f7..08daa03 100644
--- a/main.c
+++ b/main.c
@@ -1030,7 +1030,7 @@ int compress(double* data, double *m, int nrow, int ncol)
   wn = (double*) calloc(ncol * ncol, sizeof(double));
 
   get_wavelet_matrix(wn, ncol);
-  transpose(wn, nrow, ncol);
+  transpose(wn, ncol, ncol);
 
   matmul(c, nrow, ncol, wn, ncol, ncol, m);
 
diff --git a/tooling/mammoth-lagrit b/tooling/mammoth-lagrit
new file mode 100755
index 0000000..f17cfd1
--- /dev/null
+++ b/tooling/mammoth-lagrit
@@ -0,0 +1,28 @@
+#!/usr/bin/env bash
+
+echo "-------------------------------------------------------- "
+echo "mammoth lagrit: a lagrit wrapper to create exodus meshes"
+echo "-------------------------------------------------------- "
+
+sleep 5
+
+file=`ls -1 ../output/*.inp | tail -1`
+echo "[ok] using ${file} to stack mesh"
+
+cp ${file} ./layer_top.inp
+lagrit < scripts-lagrit/stack_from_top.lgi
+
+trash layer* lagrit.log lagrit.out
+
+mv -v stacked_mesh.inp scripts-lagrit/
+pushd scripts-lagrit/
+lagrit < main.lgi
+
+trash *.avs lagrit.log lagrit.out stacked_mesh.inp
+mv -v mesh.exo ../../output/
+
+popd
+
+echo "-------------------------------------------------------- "
+echo "[ok] done "
+echo "-------------------------------------------------------- "
diff --git a/tooling/scripts-lagrit/infile_get_facesets.mlgi b/tooling/scripts-lagrit/infile_get_facesets.mlgi
new file mode 100644
index 0000000..3835a85
--- /dev/null
+++ b/tooling/scripts-lagrit/infile_get_facesets.mlgi
@@ -0,0 +1,182 @@
+# ============================================================================
+# ============================================================================
+#
+# infile_get_facesets.mlgi
+# ------------------------
+#
+# Subroutine to generate user specified face sets, which are used to define
+# boundary conditions in the hydrological model ATS.
+#
+# ----------------------------------------------------------------------------
+#
+# author:
+# 	ilhan ozgen xian, iozgen@lbl.gov
+#
+# changes:
+#	February 3rd, 2020
+#
+# ============================================================================
+# ============================================================================
+
+
+
+# ============================================================================
+# -- reorder the elements --
+# ============================================================================
+
+# we prepare the mesh by rearranging the elements based on their median
+# points.
+
+cmo select CMO_PRISM
+createpts / median
+sort / CMO_PRISM / index / ascending / ikey / itetclr xmed ymed zmed
+reorder / CMO_PRISM / ikey
+cmo / DELATT / CMO_PRISM / ikey
+
+# ============================================================================
+# -- end of reorder the elements --
+# ============================================================================
+
+#
+#
+
+# ============================================================================
+# -- extract the surface mesh --
+# ============================================================================
+
+extract/surfmesh/1,0,0/cmo_tmp_surface/CMO_PRISM/external
+
+# -----------------------------------------------------------------------------
+# -- change material id of the surface mesh elements
+# -----------------------------------------------------------------------------
+cmo/select/cmo_tmp_surface
+settets/normal
+cmo/setatt/cmo_tmp_surface/itetclr 3
+
+# ============================================================================
+# -- end of extract the surface mesh: now all sides have the material id 3  --
+# ============================================================================
+
+#
+#
+
+# ============================================================================
+# -- create face sets --
+# ============================================================================
+
+# -----------------------------------------------------------------------------
+# 1. create bottom faceset
+# -----------------------------------------------------------------------------
+
+# we know that the bottom layer has the attribute layertyp -1,
+# so we select all points with this layertyp
+pset/p1/attribute/layertyp/1,0,0/-1/eq
+eltset/e1/exclusive/pset/get/p1
+cmo/setatt/cmo_tmp_surface/itetclr eltset,get,e1 1
+cmo/copy/cmo_tmp_1/cmo_tmp_surface
+cmo/DELATT/cmo_tmp_1/itetclr0
+cmo/DELATT/cmo_tmp_1/itetclr1
+cmo/DELATT/cmo_tmp_1/facecol
+cmo/DELATT/cmo_tmp_1/idface0
+cmo/DELATT/cmo_tmp_1/idelem0
+
+# from all points, select those points with layertype not equal to -1 and
+# delete them
+eltset/eall/itetclr/ge/0
+eltset/edel/not eall e1
+rmpoint/element/eltset get edel
+rmpoint/compress
+
+# -----------------------------------------------------------------------------
+# -- write the face set
+# -----------------------------------------------------------------------------
+dump/avs2/fs1_bottom.avs/cmo_tmp_1/0 0 0 2
+# -----------------------------------------------------------------------------
+
+# we will clean up after ourselves
+cmo/delete/cmo_tmp_1
+
+# -----------------------------------------------------------------------------
+# -- end of create bottom faceset
+# -----------------------------------------------------------------------------
+
+# -----------------------------------------------------------------------------
+# -- 2. create top faceset
+# -----------------------------------------------------------------------------
+
+# this time, we know that the surface layer has the attribute layertyp -2
+# and we select all points with this layertyp
+cmo/select/cmo_tmp_surface
+pset/p2/attribute/layertyp/1,0,0/-2/eq
+eltset/e2/exclusive/pset/get/p2
+cmo/setatt/cmo_tmp_surface/itetclr eltset,get,e2 2
+cmo/copy/cmo_tmp_1/cmo_tmp_surface
+cmo/DELATT/cmo_tmp_1/itetclr0
+cmo/DELATT/cmo_tmp_1/itetclr1
+cmo/DELATT/cmo_tmp_1/facecol
+cmo/DELATT/cmo_tmp_1/idface0
+cmo/DELATT/cmo_tmp_1/idelem0
+
+# from all points, select those points with layertype not equal to -2 and
+# delete them
+eltset/eall/itetclr/ge/0
+eltset/edel/not eall e2
+rmpoint/element/eltset get edel
+rmpoint/compress
+
+# -----------------------------------------------------------------------------
+# -- write the face set
+# -----------------------------------------------------------------------------
+dump/avs2/fs2_top.avs/cmo_tmp_1/0 0 0 2
+# -----------------------------------------------------------------------------
+
+# we will clean up after ourselves
+cmo/delete/cmo_tmp_1
+
+# -----------------------------------------------------------------------------
+# -- end of create top faceset
+# -----------------------------------------------------------------------------
+
+# -----------------------------------------------------------------------------
+# 3. create side faceset without directions
+# -----------------------------------------------------------------------------
+
+# material id of side sets itetclr is 3
+cmo select cmo_tmp_surface
+cmo/copy/cmo_tmp_1/cmo_tmp_surface
+
+cmo/DELATT/cmo_tmp_1/itetclr0
+cmo/DELATT/cmo_tmp_1/itetclr1
+cmo/DELATT/cmo_tmp_1/facecol
+cmo/DELATT/cmo_tmp_1/idface0
+cmo/DELATT/cmo_tmp_1/idelem0
+
+# from all points, select those points less than 3 and delete them
+eltset/edel/itetclr lt 3
+rmpoint/element/eltset get edel
+rmpoint/compress
+
+# -----------------------------------------------------------------------------
+# -- write the face set
+# -----------------------------------------------------------------------------
+dump/avs2/fs3_sides_all.avs/cmo_tmp_1/0 0 0 2
+# -----------------------------------------------------------------------------
+
+# -----------------------------------------------------------------------------
+# -- we will clean up after ourselves
+# -----------------------------------------------------------------------------
+cmo/delete/cmo_tmp_1
+# -----------------------------------------------------------------------------
+
+# -----------------------------------------------------------------------------
+# -- end of create side faceset
+# -----------------------------------------------------------------------------
+
+#
+#
+
+# ============================================================================
+# -- end of this subroutine --
+# ============================================================================
+
+finish
diff --git a/tooling/scripts-lagrit/infile_stack_from_top.mlgi b/tooling/scripts-lagrit/infile_stack_from_top.mlgi
new file mode 100644
index 0000000..b705949
--- /dev/null
+++ b/tooling/scripts-lagrit/infile_stack_from_top.mlgi
@@ -0,0 +1,13 @@
+# automatically generated stack file
+cmo/create/CMO_STACK
+cmo/select/CMO_STACK
+stack/layers/avs &
+layer1.inp 1 &
+layer_top.inp 1 / flip
+dump/avs/_tmp.inp/CMO_STACK
+# fill volumes with hex or pri
+stack/fill/CMO_PRISM/CMO_STACK
+cmo/select/CMO_PRISM
+resetpts/itp
+dump/avs/_tmp_pr.inp/CMO_PRISM
+finish
diff --git a/tooling/scripts-lagrit/main.lgi b/tooling/scripts-lagrit/main.lgi
new file mode 100644
index 0000000..795cd24
--- /dev/null
+++ b/tooling/scripts-lagrit/main.lgi
@@ -0,0 +1,93 @@
+# ============================================================================
+# ============================================================================
+#
+# Generate a stacked triangular mesh with user defined face sets for boundary
+# conditions. The surface layer mesh is extruded in the vertical direction.
+# Each subsurface layer follows the surface layer.
+#
+# ----------------------------------------------------------------------------
+#
+# author:
+# 	ilhan ozgen xian, iozgen@lbl.gov
+#
+# changes:
+#	February 3rd, 2020
+#
+# ============================================================================
+# ============================================================================
+
+#
+# ----------------------------------------------------------------------------
+#
+
+# ============================================================================
+# -- globally defined variables --
+# ============================================================================
+
+define/CMO_STACK/cmo_stack
+define/CMO_PRISM/cmo_prism
+define/STACKEDMESH/stacked_mesh.inp
+
+# ============================================================================
+# -- end of globally defined variables --
+# ============================================================================
+
+#
+#
+
+# ============================================================================
+# -- handling the top layer --
+# ============================================================================
+
+# ----------------------------------------------------------------------------
+# -- case 1: mesh the region inside this file:
+# ----------------------------------------------------------------------------
+
+# if we have no stacked mesh, LaGriT must stack several layers on top
+# of each other. in this case, uncomment the following line to call the
+# subroutine that will do this:
+# infile/infile_stack_from_top.mlgi
+
+# the result is saved in CMO_PRISM
+
+# ----------------------------------------------------------------------------
+# -- case 2: read the top layer:
+# ----------------------------------------------------------------------------
+
+# if we already have a stacked  mesh file, we can read it in, which will
+# save time. read in the file; store it in CMO_PRISM
+
+read/STACKEDMESH/CMO_PRISM
+
+# ============================================================================
+# -- end of handling the top layer
+# ============================================================================
+
+#
+#
+
+# ============================================================================
+# -- generate face sets --
+# ============================================================================
+
+infile/infile_get_facesets.mlgi
+
+# ============================================================================
+# -- end of generate face sets --
+# ============================================================================
+
+#
+#
+
+# ============================================================================
+# -- build mesh --
+# ============================================================================
+
+dump/exo/mesh.exo/CMO_PRISM///facesets &
+fs1_bottom.avs fs2_top.avs fs3_sides_all.avs
+
+# ============================================================================
+# -- end of build mesh --
+# ============================================================================
+
+finish
diff --git a/tooling/scripts-lagrit/stack_from_top.lgi b/tooling/scripts-lagrit/stack_from_top.lgi
new file mode 100644
index 0000000..eda4d12
--- /dev/null
+++ b/tooling/scripts-lagrit/stack_from_top.lgi
@@ -0,0 +1,51 @@
+define/TOP_FILE/layer_top.inp
+
+read/TOP_FILE/cmo_layer
+cmo/select/cmo_layer
+
+dump/avs/layer_top.inp cmo_layer
+
+# -- shallow soil
+
+math/sub/cmo_layer/zic/1,0,0/cmo_layer/zic/ 0.10
+dump/avs/layer1.inp/cmo_layer
+
+math/sub/cmo_layer/zic/1,0,0/cmo_layer/zic/ 0.55
+dump/avs/layer2.inp/cmo_layer
+
+# -- weathered shale
+
+math/sub/cmo_layer/zic/1,0,0/cmo_layer/zic/ 2.0
+dump/avs/layer3.inp/cmo_layer
+
+math/sub/cmo_layer/zic/1,0,0/cmo_layer/zic/ 2.35
+dump/avs/layer4.inp/cmo_layer
+
+# -- fractured shale
+
+math/sub/cmo_layer/zic/1,0,0/cmo_layer/zic/ 15.0
+dump/avs/layer5.inp/cmo_layer
+
+math/sub/cmo_layer/zic/1,0,0/cmo_layer/zic/ 20.0
+dump/avs/layer6.inp/cmo_layer
+
+cmo/create/cmo_prism
+cmo/create/cmo_stack
+cmo/select/cmo_stack
+
+stack/layers/avs/ &
+	layer6.inp 6 &
+	layer5.inp 5 &
+	layer4.inp 4 &
+	layer3.inp 3 &
+	layer2.inp 2 &
+	layer1.inp 1 &
+	layer_top.inp 1 / flip
+
+stack/fill/cmo_prism/cmo_stack
+cmo/select/cmo_prism
+resetpts/itp
+
+dump/avs/stacked_mesh.inp/cmo_prism
+
+finish