From 482421c42749b95071776a082049eb6a9916a146 Mon Sep 17 00:00:00 2001
From: Pascal FREY <>
Date: Wed, 1 Aug 2018 11:07:09 +0200
Subject: [PATCH] modified for mass computation

---
 sources/inout.c   | 26 --------------------------
 sources/nstokes.c |  2 +-
 sources/packing.c | 35 +++++++++++++++++++++++++++++++----
 3 files changed, 32 insertions(+), 31 deletions(-)

diff --git a/sources/inout.c b/sources/inout.c
index ae8e11f..05c9da2 100755
--- a/sources/inout.c
+++ b/sources/inout.c
@@ -3,11 +3,6 @@
 #include "libmesh5.h"
 
 
-static inline double mass_2d(double *a,double *b,double *c) {
-  return( ((b[0]-a[0])*(c[1]-a[1]) - (b[1]-a[1])*(c[0]-a[0])) / 2.0 );
-}
-
-
 /* read mesh */
 int loadMesh(NSst *nsst) {
   pMat       pm;
@@ -118,19 +113,6 @@ int loadMesh(NSst *nsst) {
     for (k=1; k<=nsst->info.nt; k++) {
       pt1 = &nsst->mesh.tria[k];
       GmfGetLin(inm,GmfTriangles,&pt1->v[0],&pt1->v[1],&pt1->v[2],&pt1->ref);
-      /* compute initial mass */
-      if ( nsst->info.dim == 2 && nsst->sol.nmat ) {
-        a = &nsst->mesh.point[pt1->v[0]].c[0];
-        b = &nsst->mesh.point[pt1->v[1]].c[0];
-        c = &nsst->mesh.point[pt1->v[2]].c[0];
-        for (i=0; i<nsst->sol.nmat; i++) {
-          pm = &nsst->sol.mat[i];
-          if ( pm->ref == pt1->ref ) {
-            pm->mass += mass_2d(a,b,c);
-            break;
-          }
-        }
-      }
     }
   }
   /* read mesh tetrahedra */
@@ -151,14 +133,6 @@ int loadMesh(NSst *nsst) {
     if ( nsst->info.nt )  fprintf(stdout,", %d triangles",nsst->info.nti);
     if ( nsst->info.ne )  fprintf(stdout,", %d tetrahedra",nsst->info.nei);
     fprintf(stdout,"\n");
-    if ( nsst->sol.nmat > 0 ) {
-      printf("    mass: ");
-      for (i=0; i<nsst->sol.nmat; i++) {
-        pm = &nsst->sol.mat[i];
-        printf("%g, ",pm->mass);
-      }
-      fprintf(stdout,"\n");
-    }
   }
 
   return(1);
diff --git a/sources/nstokes.c b/sources/nstokes.c
index c3568e8..908fd49 100755
--- a/sources/nstokes.c
+++ b/sources/nstokes.c
@@ -523,7 +523,7 @@ int main(int argc,char **argv) {
   }
 
   /* packing mesh if needed */
-  if ( nsst.sol.nmat ) {
+  if ( nsst.sol.nmat > 0 ) {
     ier = nsst.info.dim == 2 ? pack_2d(&nsst) : pack_3d(&nsst);
     if ( ier == 0 ) {
       if ( nsst.info.verb != '0' )  fprintf(stdout," # Packing error.\n");
diff --git a/sources/packing.c b/sources/packing.c
index 343e68d..768562e 100644
--- a/sources/packing.c
+++ b/sources/packing.c
@@ -1,6 +1,11 @@
 #include "nstokes.h"
 
 
+static inline double mass_2d(double *a,double *b,double *c) {
+  return( ((b[0]-a[0])*(c[1]-a[1]) - (b[1]-a[1])*(c[0]-a[0])) / 2.0 );
+}
+
+
 /* compactify mesh structure */
 int pack_3d(NSst *nsst) {
   pTetra    pe;
@@ -184,17 +189,39 @@ int pack_2d(NSst *nsst) {
   pTria     pt;
   pEdge     pa;
   pPoint    ppt;
-  double    nu,rho,w[2];
+  pMat      pm;
+  double   *a,*b,*c,mass,nu,rho,w[2];
   int      *prm,i,k,nf,id;
 
   /* check if compression needed */
   nf = 0;
   for (k=1; k<=nsst->info.nti; k++) {
     pt = &nsst->mesh.tria[k];
-    if ( getMat(&nsst->sol,pt->ref,&nu,&rho) ) {
-      nf++;
-      for (i=0; i<3; i++)  nsst->mesh.point[pt->v[i]].new = pt->v[i];
+    a  = &nsst->mesh.point[pt->v[0]].c[0];
+    b  = &nsst->mesh.point[pt->v[1]].c[0];
+    c  = &nsst->mesh.point[pt->v[2]].c[0];
+
+    /* compute mass */
+    for (i=0; i<nsst->sol.nmat; i++) {
+      pm = &nsst->sol.mat[i];
+      if ( pm->ref == pt->ref ) {
+        pm->mass += mass_2d(a,b,c);
+        nf++;
+        /* for renumbering purposes */
+        for (i=0; i<3; i++)  nsst->mesh.point[pt->v[i]].new = pt->v[i];
+        break;
+      }
+    }
+  }
+  if ( nsst->info.verb != '0' ) {
+    fprintf(stdout,"    mass: ");
+    mass = 0.0;
+    for (i=0; i<nsst->sol.nmat; i++) {
+      pm = &nsst->sol.mat[i];
+      fprintf(stdout," %g [%d], ",pm->mass,pm->ref);
+      mass += pm->mass;
     }
+    fprintf(stdout,"  total: %g\n",mass);
   }
   if ( nf == nsst->info.nti )  return(-1);