src/MULTIPR/MULTIPR_GaussLoc.cxx

   1 // Project MULTIPR, IOLS WP1.2.1 - EDF/CS
   2 // Partitioning/decimation module for the SALOME v3.2 platform
   3
   4 /**
   5  * \file    MULTIPR_GaussLoc.cxx
   6  *
   7  * \brief   see MULTIPR_GaussLoc.hxx
   8  *
   9  * \author  Olivier LE ROUX - CS, Virtual Reality Dpt
  10  *
  11  * \date    01/2007
  12  */
  13
  14 //*****************************************************************************
  15 // Includes section
  16 //*****************************************************************************
  17
  18 #include "MULTIPR_GaussLoc.hxx"
  19 #include "MULTIPR_Utils.hxx"
  20 #include "MULTIPR_Exceptions.hxx"
  21
  22 #include <iostream>
  23
  24 using namespace std;
  25
  26
  27 namespace multipr
  28 {
  29
  30
  31 //*****************************************************************************
  32 // Class GaussLoc implementation
  33 //*****************************************************************************
  34
  35 GaussLoc::GaussLoc()
  36 {
  37         mRefCoo    = NULL;
  38         mGaussCoo  = NULL;
  39         mWeight    = NULL;
  40
  41         reset();
  42 }
  43
  44
  45 GaussLoc::GaussLoc(const GaussLoc& pGaussLoc)
  46 {
  47         mRefCoo   = NULL;
  48         mGaussCoo = NULL;
  49         mWeight   = NULL;
  50
  51         strcpy(mName, pGaussLoc.mName);
  52
  53         mGeom     = pGaussLoc.mGeom;
  54         mDim      = pGaussLoc.mDim;
  55         mNumNodes = pGaussLoc.mNumNodes;
  56         mNumGauss = pGaussLoc.mNumGauss;
  57
  58         if (mDim != (mGeom / 100)) throw IllegalStateException("", __FILE__, __LINE__);
  59         if (mNumNodes != (mGeom % 100)) throw IllegalStateException("", __FILE__, __LINE__);
  60
  61         if (pGaussLoc.mRefCoo != NULL)
  62         {
  63                 mRefCoo = new med_float[mDim * mNumNodes];
  64                 memcpy(mRefCoo, pGaussLoc.mRefCoo, sizeof(med_float) * mDim * mNumNodes);
  65         }
  66
  67         if (pGaussLoc.mGaussCoo != NULL)
  68         {
  69                 mGaussCoo  = new med_float[mDim * mNumGauss];
  70                 memcpy(mGaussCoo, pGaussLoc.mGaussCoo, sizeof(med_float) * mDim * mNumGauss);
  71         }
  72
  73         if (pGaussLoc.mWeight != NULL)
  74         {
  75                 mWeight = new med_float[mNumGauss];
  76                 memcpy(mWeight, pGaussLoc.mWeight, sizeof(med_float) * mNumGauss);
  77         }
  78 }
  79
  80
  81 GaussLoc::~GaussLoc()
  82 {
  83         reset();
  84 }
  85
  86
  87 void GaussLoc::reset()
  88 {
  89         mName[0]  = '\0';
  90         mGeom     = MED_NONE;
  91         mDim      = 0;
  92         mNumNodes = 0;
  93         mNumGauss = 0;
  94
  95         if (mRefCoo   != NULL) { delete[] mRefCoo;   mRefCoo   = NULL; }
  96         if (mGaussCoo != NULL) { delete[] mGaussCoo; mGaussCoo = NULL; }
  97         if (mWeight   != NULL) { delete[] mWeight;   mWeight   = NULL; }
  98 }
  99
 100
 101 void GaussLoc::getCoordGaussPoints(
 102         const med_float* pCooElt,
 103         med_float* pCooGaussPoints) const
 104 {
 105         // debug
 106         //printArray2D(pCooElt, mNumNodes, mDim, "Node");
 107
 108         // WARNING: assumes TETRA10 !!!
 109         // This method is not completely generic and should be extended to support all cases.
 110         if (mGeom != MED_TETRA10) throw UnsupportedOperationException("only support TETRA10 for the moment", __FILE__, __LINE__);
 111
 112         const med_float* pt1 = pCooElt;
 113         const med_float* pt2 = pt1 + mDim;
 114         const med_float* pt3 = pt2 + mDim;
 115         const med_float* pt4 = pt3 + mDim;
 116
 117         const med_float* coeff = mGaussCoo;
 118         med_float* dest        = pCooGaussPoints;
 119
 120         // for each Gauss point
 121         for (int i = 0 ; i < mNumGauss ; i++)
 122         {
 123                 dest[0] = pt2[0] + (pt4[0] - pt2[0]) * coeff[0] + (pt1[0] - pt2[0]) * coeff[1] + (pt3[0] - pt2[0]) * coeff[2];
 124                 dest[1] = pt2[1] + (pt4[1] - pt2[1]) * coeff[0] + (pt1[1] - pt2[1]) * coeff[1] + (pt3[1] - pt2[1]) * coeff[2];
 125                 dest[2] = pt2[2] + (pt4[2] - pt2[2]) * coeff[0] + (pt1[2] - pt2[2]) * coeff[1] + (pt3[2] - pt2[2]) * coeff[2];
 126
 127                 // prepare next point
 128                 coeff += mDim;
 129                 dest += mDim;
 130         }
 131 }
 132
 133
 134 void GaussLoc::readMED(med_idt pMEDfile, med_int pIndex)
 135 {
 136         if (pMEDfile == 0) throw IOException("pMEDfile should not be NULL", __FILE__, __LINE__);
 137         if (pIndex < 1) throw IllegalArgumentException("pIndex should be >= 1", __FILE__, __LINE__);
 138
 139         reset();
 140
 141         med_err ret = MEDgaussInfo(
 142                 pMEDfile,
 143                 pIndex,
 144                 mName,
 145                 &mGeom,
 146                 &mNumGauss);
 147
 148         if (ret != 0) throw IOException("i/o error while reading Gauss localization information in MED file", __FILE__, __LINE__);
 149
 150         mDim      = mGeom / 100;
 151         mNumNodes = mGeom % 100;
 152
 153         mRefCoo   = new med_float[mDim * mNumNodes];
 154         mGaussCoo = new med_float[mDim * mNumGauss];
 155         mWeight   = new med_float[mNumGauss];
 156
 157         ret = MEDgaussLire(
 158                 pMEDfile,
 159                 mRefCoo,
 160                 mGaussCoo,
 161                 mWeight,
 162                 MED_FULL_INTERLACE,
 163                 mName);
 164
 165         if (ret != 0) throw IOException("i/o error while reading Gauss localization in MED file", __FILE__, __LINE__);
 166 }
 167
 168
 169 void GaussLoc::writeMED(med_idt pMEDfile)
 170 {
 171         if (pMEDfile == 0) throw IOException("", __FILE__, __LINE__);
 172         if (mNumGauss < 0) throw IllegalStateException("", __FILE__, __LINE__);
 173         if (mRefCoo == NULL) throw IllegalStateException("", __FILE__, __LINE__);
 174         if (mGaussCoo == NULL) throw IllegalStateException("", __FILE__, __LINE__);
 175         if (mWeight == NULL) throw IllegalStateException("", __FILE__, __LINE__);
 176         if (strlen(mName) > MED_TAILLE_NOM) throw IllegalStateException("", __FILE__, __LINE__);
 177
 178         med_err ret = MEDgaussEcr(
 179                 pMEDfile,
 180                 mGeom,
 181                 mRefCoo,
 182                 MED_FULL_INTERLACE,
 183                 mNumGauss,
 184                 mGaussCoo,
 185                 mWeight,
 186                 mName);
 187
 188         if (ret != 0) throw IOException("i/o error while writing Gauss localization", __FILE__, __LINE__);
 189 }
 190
 191
 192 ostream& operator<<(ostream& pOs, GaussLoc& pG)
 193 {
 194         pOs << "Gauss ref:" << endl;
 195         pOs << "    Name     =|" << pG.mName << "|" << endl;
 196         pOs << "    Geom     =" << pG.mGeom << endl;
 197         pOs << "    #Pt Gauss=" << pG.mNumGauss << endl;
 198
 199         pOs << "    Ref nodes coords.: (#nodes=" << pG.mNumNodes << " dim=" << pG.mDim << ")" << endl;
 200         for (int itNode = 0 ; itNode < pG.mNumNodes ; itNode++)
 201         {
 202                 pOs << "        Node " << (itNode + 1) << ": ";
 203                 for (int itDim = 0; itDim < pG.mDim ; itDim++)
 204                 {
 205                         pOs << pG.mRefCoo[itNode * pG.mDim + itDim] << " ";
 206                 }
 207                 pOs << endl;
 208         }
 209
 210         pOs << "    Gauss coords. and weight:" << endl;
 211         for (int itGauss = 0 ; itGauss < pG.mNumGauss ; itGauss++)
 212         {
 213                 pOs << "        Pt " << (itGauss+1) << ": ";
 214                 for (int itDim = 0; itDim < pG.mDim ; itDim++)
 215                 {
 216                         pOs << pG.mGaussCoo[itGauss * pG.mDim + itDim] << " ";
 217                 }
 218                 pOs << "weight=" << pG.mWeight[itGauss];
 219                 pOs << endl;
 220         }
 221
 222         /*
 223         // debug
 224         med_float res[15];
 225         pG.getCoordGaussPoints(pG.mRefCoo, res);
 226         printArray2D(res, pG.mNumGauss, pG.mDim, "Gauss pt");
 227         */
 228
 229         return pOs;
 230 }
 231
 232
 233 } // namespace multipr
 234
 235 // EOF