matrix.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:                          ###  ###  #####   ####  ##  ##
//                       ##  ##   ########  ##     ##  ## ##  ##
//                       ##  ##   ## ## ##  ###    ##     ######
//                       ##  ##   ##    ##  ##     ##  ## ##  ##
//                       #####    ##    ##  #####   ####  ##  ##
//                       ##
//
// name:           matrix_record.cpp
// author(s):      Jan Novak
// last edit:      04. 02. 2010
// language:       C, C++
// license:        This program is free software; you can redistribute it and/or modify
//                 it under the terms of the GNU Lesser 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 Lesser General Public License for more details.
//
//                 You should have received a copy of the GNU Lesser General Public License
//                 along with this program; if not, write to the Free Software
//                 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
//********************************************************************************************************

//included libraries
#include <stdio.h>
#include "matrix.h"
#include "stiffness.h"
#include "esuf.h"
#include "matrixOperations.h"
#include "problem.h"
#include "librw.h"


namespace mumech {


MatrixRecord :: MatrixRecord (const Problem *p)
{
  P = p;
  ym = pr = -1.0;
  origin = NULL;
  nRS = 0;
  Remote_strain = NULL;
  
  globHomog_displc = globHomog_strain = globHomog_stress = NULL;

  C = new double[P->give_ISO_C_RANGE()];
}

MatrixRecord :: ~MatrixRecord ()
{
  delete [] origin;
  DeleteArray2D (Remote_strain, P->give_nLC());
  DeleteArray2D (globHomog_displc, P->give_nLC());
  DeleteArray2D (globHomog_strain, P->give_nLC());
  DeleteArray2D (globHomog_stress, P->give_nLC());

  if (C)    delete [] C;
}

void MatrixRecord :: set_nlc (int n)
{
  if (nRS) _errorr("nlc set already");
  nRS = n;
  if (nRS < 1  ||  nRS > 99)  _errorr("Number of load cases should be > 0 and <100");
}

void MatrixRecord :: set_E_nu (double e, double n)
{
  if (P->is_converted_to_equivalent())
    _errorr("You can not set matrix recor when data already converted to equivalent problem.");
  
  if (this->ym >= 0.0)  _errorr("Matrix record: E is already set");
  if (this->pr >= 0.0)  _errorr("Matrix record: nu is already set");
  
  ym = e;
  pr = n;
  
  // evaluation of stifness matrix
  Stiffness::giveReducedIsoStiffMatrix (this->C, this->ym, this->pr, P->is_twodim());
}

void MatrixRecord :: set_origin (double x, double y)           {  P->check_dim(2); origin = new double[2]; origin[0] = x; origin[1] = y;                }
void MatrixRecord :: set_origin (double x, double y, double z) {  P->check_dim(3); origin = new double[3]; origin[0] = x; origin[1] = y; origin[2] = z; }

void MatrixRecord :: convert_RS_ROWtoFEEP (int lc)
{
  if (P->is_twodim()) MatrixOperations::convert2Dtensor_4ROWto3FEEP_notation (Remote_strain[lc]);
  else                MatrixOperations::convert3Dtensor_9ROWto6FEEP_notation (Remote_strain[lc]);
}

void MatrixRecord :: set_Remote_strain_for_lc (int lc, const double r1[])
{
  if (this->nRS == 0)  _errorr("Number of remote strains was not set.");
  if (this->nRS <= lc) _errorr("Remote strain ID is greather then the set number of remote strains");
  if (Remote_strain == NULL)  this->allocate_nlc_fields ();
  CopyVector (r1, Remote_strain[lc], P->give_TENS_RANGE());
  
  this->convert_RS_ROWtoFEEP (lc);
}
bool MatrixRecord :: scan_Remote_strains (FILE *stream)
{
  if (Remote_strain == NULL)  this->allocate_nlc_fields ();
  
  int l = P->give_TENS_RANGE();
  
  for (int i=0; i<P->give_nLC(); i++) {
    if (FP_scan_array (stream, l, Remote_strain[i]) != true)
      return false;
    
    this->convert_RS_ROWtoFEEP (i);
  }
  
  return true;
}
bool MatrixRecord :: scan_Remote_strains (const char *stream)
{
  if (Remote_strain == NULL)  this->allocate_nlc_fields ();
  
  int l = P->give_TENS_RANGE();
  
  for (int i=0; i<P->give_nLC(); i++) {
    if (SP_scan_array (stream, l, Remote_strain[i]) != true)
      return false;
    
    this->convert_RS_ROWtoFEEP (i);
  }
  
  return true;
}

void MatrixRecord :: set_unit_remote_strain (void)
{
  if (Remote_strain) {
    if (this->remote_strain_is_unit() == false)
      _errorr("Already given Remote strain is not equal to \"unit set\".");

  }
  else {
    // Unit remote strain se pouziva pro homogenizaci RegGrid, a musi byt proto na vsech pozicich 1.0.
    // Ikdyz tam dam 0.5, tak se stejne nic nedeje, pac to tam stejne nasobim inverzi strainu,
    // takze vlastne muze byt unit strain libovolnej, ale jen diky tomu nasobeni.
    double smyk = 1.0;

    this->set_nlc(P->give_VM_TENS_RANGE());
    this->allocate_nlc_fields();
    
    CleanArray2d (Remote_strain, nRS, P->give_TENS_RANGE());
    
    if (P->twodim) {
      Remote_strain[0][0] = 1.0;
      Remote_strain[1][1] = 1.0;
      Remote_strain[2][2] = smyk;
    }
    else {
      Remote_strain[0][0] = 1.0;
      Remote_strain[1][1] = 1.0;
      Remote_strain[2][2] = 1.0;
      Remote_strain[4][4] = smyk;
      Remote_strain[5][5] = smyk;
      Remote_strain[3][3] = smyk;
    }
  }
}

bool MatrixRecord :: remote_strain_is_unit (void) const
{
  if (this->nRS != P->give_VM_TENS_RANGE())  return false;

  double smyk = 1.0; // same as upper

  double s = GiveSumMembersArray2D ((const double **)Remote_strain, nRS, P->give_TENS_RANGE());
  if (s != (P->twodim ? 2+smyk : 3+3*smyk))  return false;

  if (P->twodim) {
    if (Remote_strain[0][0] != 1.0 ||
        Remote_strain[1][1] != 1.0 ||
        Remote_strain[2][2] != smyk   )  return false;
  }
  else {
    if (Remote_strain[0][0] != 1.0  ||
        Remote_strain[1][1] != 1.0  ||
        Remote_strain[2][2] != 1.0  ||
        Remote_strain[3][3] != smyk ||
        Remote_strain[4][4] != smyk ||
        Remote_strain[5][5] != smyk    )  return false;
  }
  return true;
}

void MatrixRecord :: giveReducedStiffMatrix (double *answer) const
{
  CopyVector(this->C, answer, P->give_ISO_C_RANGE());
}
void MatrixRecord :: giveFullStiffMatrix (double *answer) const
{
  if (P->twodim)  MatrixOperations::copy2DeshelbyTensor_reduced2full(this->C, answer);
  else            MatrixOperations::copy3DeshelbyTensor_reduced2full(this->C, answer);
}


void MatrixRecord :: allocate_nlc_fields (void)
{
  AllocateArray2D (Remote_strain,    nRS, P->give_TENS_RANGE());
  AllocateArray2D (globHomog_displc, nRS, P->give_TENS_RANGE());
}

const double* MatrixRecord :: give_globHomog_Displc (int lc, const double *coords)
{
  MatrixOperations::giveIsoTensAndVectorProduct(Remote_strain[lc], coords, globHomog_displc[lc], P->give_dimension());
  
  return globHomog_displc[lc];
}
void MatrixRecord :: give_globHomog_Displc (double *displ, int lc, const double *coords)
{
  CopyVector(this->give_globHomog_Displc(lc, coords), displ, P->give_VECT_RANGE());
}

const double* MatrixRecord :: give_globHomog_Strain (int lc)
{
  return Remote_strain[lc];
}
void MatrixRecord :: give_globHomog_Strain (double *strain, int lc)
{
  CopyVector(this->give_globHomog_Strain(lc), strain, P->give_VM_TENS_RANGE());
}

const double* MatrixRecord :: give_globHomog_Stress (int lc)
{
  if (globHomog_stress == NULL) {
    globHomog_stress = new double*[P->give_nLC()];
    memset (globHomog_stress, 0, P->give_nLC() * sizeof(double*));
  }
  
  if (globHomog_stress[lc] == NULL) {
    globHomog_stress[lc] = new double[P->give_VM_TENS_RANGE()];
    
    if (P->twodim) MatrixOperations::giveTVproduct_3is3x3to5and3  (globHomog_stress[lc], this->C, Remote_strain[lc]);
    else           MatrixOperations::giveTVproduct_6is6x6to12and6 (globHomog_stress[lc], this->C, Remote_strain[lc]);
  }
  
  return globHomog_stress[lc];
}
void MatrixRecord :: give_globHomog_Stress (double *stress, int lc)
{
  CopyVector(this->give_globHomog_Stress(lc), stress, P->give_VM_TENS_RANGE());
}







void MatrixRecord :: printYourself (const char *notice)
{
  if (notice)  printf( "%s", notice );
  printf( "Young's modulus: %e\n", ym );
  printf( "Poisson's ratio: %e\n", pr );
  //eshSolu.printIsoTensor( C, "Infinite medium stifness matrix:\n", _STANDARD_SYNTAX_ );
  
  return ;
} // end of function: printMtrxRecord( )


// ********************************************************************************************************
// **************************   POINT   *******************************************************************
// ********************************************************************************************************

Point :: Point ()
{
  //position = PPF_VOID;
  //this->set_nlc(p, nlc);
  
  displacement = strain = stress = NULL;
}

Point :: ~Point ()
{
  DeleteArray2D (displacement,  nlc);
  DeleteArray2D (strain,        nlc);
  DeleteArray2D (stress,        nlc);
}


void Point :: set_nlc (const Problem *p, int nlc)
{
  Point::nlc = nlc;
  
  AllocateArray2D (displacement, nlc, p->give_VECT_RANGE());
  AllocateArray2D (strain,       nlc, p->give_VM_TENS_RANGE());
  AllocateArray2D (stress,       nlc, p->give_VM_TENS_RANGE());
}



} // end of namespace mumech

/*end of file*/