[cig-commits] r13819 - cs/cigma/trunk/src

[luis at geodynamics.org]

Author: luis
Date: 2009-01-12 14:09:39 -0800 (Mon, 12 Jan 2009)
New Revision: 13819

Modified:
   cs/cigma/trunk/src/Residuals.cpp
   cs/cigma/trunk/src/Residuals.h
Log:
Updated Residuals class to compute max cell diameter (moved here from mesh_info.cc)

Modified: cs/cigma/trunk/src/Residuals.cpp
===================================================================
--- cs/cigma/trunk/src/Residuals.cpp	2009-01-12 22:09:38 UTC (rev 13818)
+++ cs/cigma/trunk/src/Residuals.cpp	2009-01-12 22:09:39 UTC (rev 13819)
@@ -1,15 +1,13 @@
+#include "Residuals.h"
 #include <iostream>
 #include <cassert>
 #include <cmath>
-#include "Residuals.h"
 
 using namespace std;
 using namespace cigma;
 using boost::shared_ptr;
 
 
-// ----------------------------------------------------------------------------
-
 Residuals::Residuals()
 {
     nel = 0;
@@ -17,8 +15,9 @@
     volume = 0;
 
     global_error = 0;
+    total_volume = 0;
     max_residual = 0;
-    total_volume = 0;
+    max_diameter = 0;
 }
 
 Residuals::~Residuals()
@@ -42,7 +41,7 @@
 
     if (meshPart)
     {
-        nel = meshPart->connect->n_cells();
+        nel = meshPart->n_cells();
 
         if (epsilon != 0) { delete [] epsilon; }
         if (volume != 0) { delete [] volume; }
@@ -92,19 +91,120 @@
     global_error += cell_residual * cell_residual;
 }
 
-double Residuals::max()
+void Residuals::computeMaxResidual()
 {
+    assert(meshPart);
+    assert(nel > 0);
+
+    double maxres = 0.0;
+    for (int e = 0; e < nel; e++)
+    {
+        double res = epsilon[e];
+        if (res > maxres)
+        {
+            maxres = res;
+        }
+    }
+    max_residual = maxres;
+}
+
+static double dist(double dx, double dy, double dz)
+{
+    return sqrt(dx*dx + dy*dy + dz*dz);
+}
+
+static double BoxDiameter(double bounds[6])
+{
+    const double x0 = bounds[0];
+    const double x1 = bounds[1];
+    const double y0 = bounds[2];
+    const double y1 = bounds[3];
+    const double z0 = bounds[4];
+    const double z1 = bounds[5];
+
+    assert(x0 <= x1);
+    assert(y0 <= y1);
+    assert(z0 <= z1);
+
+    return dist(x1-x0, y1-y0, z1-z0);
+}
+
+void Residuals::computeMaxCellDiameter()
+{
+    assert(meshPart);
+    assert(meshPart->coords);
+    assert(meshPart->connect);
+    assert(meshPart->getCellType() != Cell::NONE);
+    assert(nel > 0);
+
+    double bounds[6];
+    int npts, ncells;
+    int e, emax;
+    double h, hmax;
+
+    shared_ptr<Cell> cell = Cell::New(meshPart->getCellType());
+    hmax = 0.0;
+    emax = -1;
+
+    for (e = 0; e < nel; e++)
+    {
+        meshPart->getCell(e, *cell);
+        cell->getBounds(bounds);
+        h = BoxDiameter(bounds);
+        if (h > hmax)
+        {
+            hmax = h;
+            emax = e;
+        }
+    }
+    max_diameter = hmax;
+}
+
+void Residuals::computeVolume()
+{
+    assert(meshPart);
+    assert(meshPart->coords);
+    assert(meshPart->connect);
+    assert(meshPart->getCellType() != Cell::NONE);
+
+    shared_ptr<Cell> cell = Cell::New(meshPart->getCellType());
+
+    double total = 0.0;
+    for (int e = 0; e < nel; e++)
+    {
+        double vol = cell->volume();
+        volume[e] = vol;
+        total += vol;
+    }
+    total_volume = total;
+}
+
+double Residuals::norm_Linf() const
+{
     return max_residual;
 }
 
-double Residuals::L2()
+double Residuals::norm_L2() const
 {
-    return sqrt(global_error);  // XXX: generalizes to norm_pow(global_error)
+    // XXX: generalizes to norm_pow(global_error)
+    return sqrt(global_error);
 }
 
-double Residuals::relative_squared_L2()
+double Residuals::relative_squared_norm_L2() const
 {
     return global_error/total_volume;
 }
 
-// ----------------------------------------------------------------------------
+double Residuals::relative_norm_L2() const
+{
+    return sqrt(global_error/total_volume);
+}
+
+double Residuals::max_cell_diameter()
+{
+    if (max_diameter <= 0)
+    {
+        this->computeMaxCellDiameter();
+    }
+    return max_diameter;
+}

Modified: cs/cigma/trunk/src/Residuals.h
===================================================================
--- cs/cigma/trunk/src/Residuals.h	2009-01-12 22:09:38 UTC (rev 13818)
+++ cs/cigma/trunk/src/Residuals.h	2009-01-12 22:09:39 UTC (rev 13819)
@@ -22,11 +22,15 @@
     void resetAccumulator();
     void update(int e, double cell_residual, double cell_volume);
 
-    double max();
-    double L2();
-    double relative_squared_L2();
+    void computeMaxResidual();
+    void computeMaxCellDiameter();
+    void computeVolume();
 
-    //void write(const char *filename); // XXX
+    double norm_Linf() const;
+    double norm_L2() const;
+    double relative_squared_norm_L2() const;
+    double relative_norm_L2() const;
+    double max_cell_diameter();
 
 public:
 
@@ -36,10 +40,10 @@
     double *epsilon;
     double *volume;
 
-
     double global_error;
+    double total_volume;
     double max_residual;
-    double total_volume;
+    double max_diameter;
 
 };