femmisc.cpp

FreeFEM is an implementation of the GFEM language dedicated to the finite element method. It provid

// Emacs will be in -*- Mode: c++ -*-
//
// ********** DO NOT REMOVE THIS BANNER **********
//
// SUMMARY: Language for a Finite Element Method
//
// AUTHORS:  C. Prud'homme
// ORG    :          
// E-MAIL :  [e-mail]
//
// ORIG-DATE:     June-94
// LAST-MOD: 24-Oct-01 at 18:43:09 by Christophe Prud'homme
//
// DESCRIPTION: 
/*
  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2 of the License, or
  (at your option) any later version.
  
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
  
  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

*/
// DESCRIP-END.
//

#ifdef __GNUG__
#pragma implementation
#endif

#include <cassert>

#include <femGraphic.hpp>
#include <femMisc.hpp>

namespace fem
{
creal sqrtofminus1 (0.F, 1.F);
int complexmode = 1;
//creal sqrtofminus1 = 0; int complexmode=0;
// Complex calculation: change this   /* C2R */

void 
myassert (int what)
{
  if (!what)
    assert (what);
}

float 
norm2 (const float &a)
{
  return a > 0 ? a : -a;
}
float 
imagpart (const float &a)
{
  return a;
}
float 
realpart (const float &a)
{
  return a;
}

std::ostream & operator << (std::ostream & os, const Complex & a)
{
  os << a.re << " " << a.im << "    ";
  return os;
}
std::istream & operator >> (std::istream & os, Complex & a)
{
  os >> a.re >> a.im;
  return os;
}

std::ostream & operator << (std::ostream & os, cvect & a)
{
  for (int i = 0; i < N; i++)
    os << a[i] << "  ";
  return os;
}
std::ostream & operator << (std::ostream & os, cmat & a)
{
  for (int i = 0; i < N; i++)
    for (int j = 0; j < N; j++)
      os << a (i, j) << "  ";
  return os;
}

Complex 
id(const Complex& x)
{
  Complex c (0.F);
  if ((x.im != 0.F)||(x.re != 0.F))
    c.re = 1.F;
  return c;
}
cmat 
id (const cvect& a)
{
  cmat c (0.F);
  for (int i = 0; i < N; i++)
    if (norm2 (a.val[i]) != 0.F)
      c (i, i) = 1.F;
  return c;
}

void 
cgauss (cmat & a, cvect & x)
{
  int i, j, k;
  ccreal s, s1;
  float smin = (float)1.0e9, eps = (float)1.0e-9;
  int lu = 1;     /* !! */

  if (lu)
    for (i = 0; i < N; i++)
      {
  for (j = 0; j <= i; j++)
    {
      s = 0.F;
      for (k = 0; k < j; k++)
        s += a (i, k) * a (k, j);
      a (i, j) -= s;

    }
  for (j = i + 1; j < N; j++)
    {
      s = 0.F;
      for (k = 0; k < i; k++)
        s += a (i, k) * a (k, j);
      s1 = a (i, i);
      if (norm2 (s1) < smin)
        smin = norm2 (s1);
      if (norm2 (s1) < eps)
        s1 = eps;
      a (i, j) = (a (i, j) - s) / s1;
    }
      }
  for (i = 0; i < N; i++)
    {
      s = 0.F;
      for (k = 0; k < i; k++)
  s += a (i, k) * x[k];
      x[i] = (x[i] - s) / a (i, i);
    }
  for (i = N - 1; i >= 0; i--)
    {
      s = 0.F;
      for (k = i + 1; k < N; k++)
  s += a (i, k) * x[k];
      x[i] -= s;
    }
//      a.lu = 1;
  //      return smin;    
}

void
Acreal::init (long newSize)
{
  myassert (!(szz || cc));
  szz = newSize;
  cc = new creal[szz];
  if (!cc)
    erreur ("Out of Memory");
  for (int i = 0; i < szz; i++)
    cc[i] = 0;
}
Acreal::Acreal (long sz)
{
  cc = 0;
  if (sz > 0)
    {
      cc = new creal[sz];
      if (!cc)
  erreur ("Out of Memory");
    }
  for (int i = 0; i < sz; i++)
    cc[i] = 0;
  szz = sz;
}
Acreal::Acreal (const Acreal & a)
{
  cc = 0;
  if (a.szz > 0)
    {
      szz = a.szz;
      cc = new creal[szz];
      if (!cc)
  erreur ("Out of Memory");
      else
  for (int i = 0; i < szz; i++)
    cc[i] = a.cc[i];
    }
  else
    {
      cc = NULL;
      szz = 0;
    }
}
Aint::Aint (long sz)
{
  cc = 0;
  if (sz > 0)
    {
      cc = new int[sz];
      if (!cc)
  erreur ("Out of Memory");
      for (int i = 0; i < sz; i++)
  cc[i] = 0;
    }
  szz = sz;
}
Aint::Aint (const Aint & a)
{
  cc = 0;
  if (a.szz > 0)
    {
      szz = a.szz;
      cc = new int[szz];
      if (!cc)
  erreur ("Out of Memory");
      else
  for (int i = 0; i < szz; i++)
    cc[i] = a.cc[i];
    }
  else
    {
      cc = NULL;
      szz = 0;
    }
}
void
Aint::init (long newSize)
{
  myassert (!(szz || cc));
  szz = newSize;
  cc = new int[szz];
  if (!cc)
    erreur ("Out of Memory");
  for (int i = 0; i < szz; i++)
    cc[i] = 0;
}
Acvect::Acvect (long sz)
{
  cc = 0;
  if (sz > 0)
    {
      cc = new cvect[sz];
      if (!cc)
  erreur ("Out of Memory");
      for (int i = 0; i < sz; i++)
  cc[i] = 0.F;
    }
  szz = sz;
}

Acvect::Acvect (const Acvect & a)
{
  cc = 0;
  if (a.szz > 0)
    {
      szz = a.szz;
      cc = new cvect[szz];
      if (!cc)
  erreur ("Out of Memory");
      else
  for (int i = 0; i < szz; i++)
    cc[i] = a.cc[i];
    }
  else
    {
      cc = NULL;
      szz = 0;
    }
}
void
Acvect::init (long newSize)
{
  myassert (!(szz || cc));
  szz = newSize;
  cc = new cvect[szz];
  if (!cc)
    erreur ("Out of Memory");
  else
    for (int i = 0; i < szz; i++)
      cc[i] = 0.F;
}
Acmat::Acmat (long sz)
{
  cc = 0;
  if (sz > 0)
    {
      cc = new cmat[sz];
      if (!cc)
  erreur ("Out of Memory");
      for (int i = 0; i < sz; i++)
  cc[i] = 0.F;
    }
  szz = sz;
}
Acmat::Acmat (const Acmat & a)
{
  cc = 0;
  if (a.szz > 0)
    {
      szz = a.szz;
      cc = new cmat[szz];
      if (!cc)
  erreur ("Out of Memory");
      else
  for (int i = 0; i < szz; i++)
    cc[i] = a.cc[i];
    }
  else
    {
      cc = NULL;
      szz = 0;
    }
}
void
Acmat::init (long newSize)
{
  myassert (!(szz || cc));
  szz = newSize;
  cc = new cmat[szz];
  if (!cc)
    erreur ("Out of Memory");
  else
    for (int i = 0; i < szz; i++)
      cc[i] = 0.F;
}
AAcmat::AAcmat (long sz)
{
  cc = 0;
  if (sz > 0)
    {
      cc = new Acmat[sz];
      if (!cc)
  erreur ("Out of Memory");
      for (int i = 0; i < sz; i++)
  cc[i] = 0;
    }
  szz = sz;
}
AAcmat::AAcmat (const AAcmat & a)
{
  cc = 0;
  if (a.szz > 0)
    {
      szz = a.szz;
      cc = new Acmat[szz];
      if (!cc)
  erreur ("Out of Memory");
      else
  for (int i = 0; i < szz; i++)
    cc[i] = a.cc[i];
    }
  else
    {
      cc = NULL;
      szz = 0;
    }
}
void
AAcmat::init (long newSize)
{
  myassert (!(szz || cc));
  szz = newSize;
  cc = new Acmat[szz];
  if (!cc)
    erreur ("Out of Memory");
  for (int i = 0; i < szz; i++)
    cc[i] = 0;
}
AAcreal::AAcreal (long sz)
{
  cc = 0;
  if (sz > 0)
    {
      cc = new Acreal[sz];
      if (!cc)
  erreur ("Out of Memory");
      for (int i = 0; i < sz; i++)
  cc[i] = 0;
    };
  szz = sz;
}
AAcreal::AAcreal (const AAcreal & a)
{
  cc = 0;
  if (a.szz > 0)
    {
      szz = a.szz;
      cc = new Acreal[szz];
      if (!cc)
  erreur ("Out of Memory");
      else
  for (int i = 0; i < szz; i++)
    cc[i] = a.cc[i];
    }
  else
    {
      cc = NULL;
      szz = 0;
    }
}
void
AAcreal::init (long newSize)
{
  myassert (!(szz || cc));
  szz = newSize;
  cc = new Acreal[szz];
  if (!cc)
    erreur ("Out of Memory");
  for (int i = 0; i < szz; i++)
    cc[i] = 0;
}
}