melement.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:                          ###  ###  #####   ####  ##  ##
//                       ##  ##   ########  ##     ##  ## ##  ##
//                       ##  ##   ## ## ##  ###    ##     ######
//                       ##  ##   ##    ##  ##     ##  ## ##  ##
//                       #####    ##    ##  #####   ####  ##  ##
//                       ##
//
// name:           element.cpp
// author(s):      Ladislav Svoboda
// 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.
//********************************************************************************************************

#include "melement.h"

#include "problem.h"
#include "mesh.h"
#include "mnode.h"
#include "matrixOperations.h"

#include "arrays.h"


namespace mumech {

mElement :: mElement (long i, const Mesh *m)
{
  id = i;
  M = m;
  
  nnodes = 0;
  nodes = NULL;
  type = 0;
  region = -3;
  volume = -1;
  
  displc = strain = stress = NULL;
  energy = NULL;
  eqstress = NULL;

  error = NULL;
}
mElement :: mElement (long i, const Mesh *m, const mElement *src)
{
  id = i;
  M = m;
  
  type = 0;
  this->set_type(src->type);
  CopyVector(src->nodes, this->nodes, this->nnodes);
  
  region = -3;
  volume = -1;

  displc = strain = stress = NULL;
  if (src->displc != NULL  ||  src->strain != NULL  ||  src->stress != NULL)  _errorr("results data not supported in copy constructor");
  energy = NULL;
  if (src->energy != NULL)  _errorr("results data not supported in copy constructor");
  eqstress = NULL;
  if (src->eqstress != NULL)  _errorr("results data not supported in copy constructor");

  error = NULL;
  if (src->error != NULL)  _errorr("results data not supported in copy constructor");

}


mElement :: ~mElement ()
{
  delete [] nodes;
  
  int nP  = M->npa;
  int nLC = M->nLC;
  
  DeleteArray3D (displc, nP, nLC);
  DeleteArray3D (strain, nP, nLC);
  DeleteArray3D (stress, nP, nLC);
  
  DeleteArray2D (energy,   nP);
  DeleteArray2D (eqstress, nP);
  DeleteArray2D (error,    nP);
}

void mElement :: set_type (long t)
{
  if (type != 0)  _errorr("");
  type = t;
  
  if      (type ==  5)  nnodes = 3;
  else if (type ==  9)  nnodes = 4;
  else if (type == 10)  nnodes = 4;
  else if (type == 12)  nnodes = 8;
  else _errorr("");
  
  nodes = new long[nnodes];
}

void mElement :: compute_centroids (void)
{
  if (nodes==NULL)  errol;
  
  if (M->coll) { //2D
    if (nnodes == 3) {
      for (int i=0; i<3; i++)
    centroids[i] = ( M->Nodes[nodes[0]]->coords[i] +
             M->Nodes[nodes[1]]->coords[i] +
             M->Nodes[nodes[2]]->coords[i]   ) / 3.0;
    }
    else if (nnodes == 4) {
      for (int i=0; i<3; i++)
    centroids[i] = ( M->Nodes[nodes[0]]->coords[i] +
             M->Nodes[nodes[1]]->coords[i] +
             M->Nodes[nodes[2]]->coords[i] +
             M->Nodes[nodes[3]]->coords[i]   ) / 4.0;
    }
    else errol;
  }
  else { // 3D
    if (nnodes == 8) {
      for (int i=0; i<3; i++)
    centroids[i] = ( M->Nodes[nodes[0]]->coords[i] + M->Nodes[nodes[4]]->coords[i] + 
             M->Nodes[nodes[1]]->coords[i] + M->Nodes[nodes[5]]->coords[i] + 
             M->Nodes[nodes[2]]->coords[i] + M->Nodes[nodes[6]]->coords[i] + 
             M->Nodes[nodes[3]]->coords[i] + M->Nodes[nodes[7]]->coords[i]   ) / 8.0;
    }
    else errol;
  }
}


void mElement :: allocate_fields (int pid)
{
  int nLC = M->nLC;
  
  if (M->el_dspl_rf) { if (!displc)  displc = Allocate1Ddpp (M->npa);  if (M->el_dspl_rf[pid]) { if (!displc[pid])  displc[pid] = Allocate1Ddp(nLC);  for (int i=0; i<nLC; i++) { if (M->el_dspl_rf[pid][i] && displc[pid][i] == NULL)  displc[pid][i] = new double[M->P[pid]->ndisplc_one()]; }  }  }
  if (M->el_strn_rf) { if (!strain)  strain = Allocate1Ddpp (M->npa);  if (M->el_strn_rf[pid]) { if (!strain[pid])  strain[pid] = Allocate1Ddp(nLC);  for (int i=0; i<nLC; i++) { if (M->el_strn_rf[pid][i] && strain[pid][i] == NULL)  strain[pid][i] = new double[M->P[pid]->nstrain_one()]; }  }  }
  if (M->el_strs_rf) { if (!stress)  stress = Allocate1Ddpp (M->npa);  if (M->el_strs_rf[pid]) { if (!stress[pid])  stress[pid] = Allocate1Ddp(nLC);  for (int i=0; i<nLC; i++) { if (M->el_strs_rf[pid][i] && stress[pid][i] == NULL)  stress[pid][i] = new double[M->P[pid]->nstrain_one()]; }  }  }
}

double mElement :: give_volume(void)
{
  if (nodes==NULL)  errol;
  
  volume = M->give_regular_element_volume();
  
  // compute volume
  if (volume <= 0.0) {
    if (M->coll) { //2D
      if (nnodes!=3) errol;
      VectoR normal;
      
      normal.beVectProduct (&M->Nodes[nodes[0]]->coords,
                &M->Nodes[nodes[1]]->coords,
                &M->Nodes[nodes[2]]->coords);
      
      volume = 0.5 * normal.give_length();
    }
    else { // 3D
      if (nnodes!=3) errol;
      errol;
    }
  }
  
  return volume;
}

double mElement :: give_energy (int pid, int lc)
{
  if (energy      == NULL)  energy      = Allocate1Ddp(M->npa);
  if (energy[pid] == NULL)  energy[pid] = Allocate1Ddz(M->nLC);
  
  // compute energy
  if (energy[pid][lc] == 0.0) {
    
    bool twodim = M->give_twodim();
    
    // mesh and problem dimension is same
    if ((M->coll == 0  &&  M->P[pid]->give_dimension() == 3) ||
        (M->coll == 3  &&  M->P[pid]->give_dimension() == 2)) {
      
      if (twodim)  energy[pid][lc] = MatrixOperations::giveTTproduct_1is2x2to3and2x2to3_SS(stress[pid][lc], strain[pid][lc]);
      else errol;
    }
    // problem is 3d, mesh is 2d => planestrain is supposed
    else if ((M->coll == 3  &&  M->P[0]->give_dimension() == 3)) {
      errol;
      //   double strain2d[3], stress2d[3];
      //   copy3Dto2Dtensors_FEEPreduced (strain[lc], strain2d);
      //   copy3Dto2Dtensors_FEEPreduced (stress[lc], stress2d);
      // 
      //   energy[lc] = give_VectorVectorProduct (strain2d, stress2d, 3);
    }
    else errol;
    
    //
    energy[pid][lc] *= this->give_volume();
  }
  
  return energy[pid][lc];
}


void mElement :: save_error (double val, int pid, int lc)
{
  if (error      == NULL)  error      = Allocate1Ddp (M->npa);
  if (error[pid] == NULL)  error[pid] = Allocate1Ddz (M->nLC);
  
  error[pid][lc] = val;
}

void mElement :: giveReducedStiffMatrix (double *C, int pid) const
{
  const Problem *prob = dynamic_cast <const Problem*>(M->P[pid]);
  if (prob == NULL) errol;
  
  if      (this->region <  -1) errol;
  else if (this->region == -1) prob->matrix_giveReducedStiffMatrix(C);
  else                         prob->give_inclusion(region)->give_StiffnessMatrixReduced(C);
}


void mElement :: give_strain (double *strain, int pid, int lc) const
{
  if (this->strain == NULL || this->strain[pid] == NULL || this->strain[pid][lc] == NULL)  errol;
  if (strain == NULL)  errol;
  
  // mesh and problem dimension is same
  if ( M->equal_dimensions(pid) ) {
    CopyVector (this->strain[pid][lc], strain, M->P[pid]->nstrain_one());
  }
  //        mesh is 2d  and    problem is 3d                         => planestrain is supposed
  else if ((M->coll == 3  &&  M->P[pid]->give_dimension() == 3)) {
    MatrixOperations::copy3Dto2Dtensors_FEEPreduced (this->strain[pid][lc], strain);
  }
  else errol;
}
void mElement :: give_stress (double *stress, int pid, int lc) const
{
  if (this->stress == NULL || this->stress[pid] == NULL || this->stress[pid][lc] == NULL)  errol;
  if (stress == NULL)  errol;

  // mesh and problem dimension is same
  if ( M->equal_dimensions(pid) ) {
    CopyVector(this->stress[pid][lc], stress, M->P[pid]->nstrain_one());
  }
  //        mesh is 2d  and    problem is 3d                         => planestrain is supposed
  else if ((M->coll == 3  &&  M->P[pid]->give_dimension() == 3)) {
    MatrixOperations::copy3Dto2Dtensors_FEEPreduced (this->stress[pid][lc], stress);
  }
  else errol;
}


// ///
// double mElement :: give_equiv_stress (int lc)
// {
//   if (eqstress == NULL) {
//     int nLC = M->P[0]->give_nLC();
//     eqstress = new double[nLC];
//     memset (eqstress, 0, nLC*sizeof(double));
//   }
// 
//   // compute eqstress
//   if (eqstress[lc] == 0.0) {
//     // mesh and problem dimension is same
//     if ((M->coll == 0  &&  M->P[0]->give_dimension() == 3) ||
//         (M->coll == 3  &&  M->P[0]->give_dimension() == 2)) {
//         errol;
//       //eqstress[lc] = give_VectorVectorProduct (strain[lc], stress[lc], M->P->nstrain());
//     }
//     // problem is 3d, mesh is 2d => planestrain is supposed
//     else if ((M->coll == 3  &&  M->P[0]->give_dimension() == 3)) {
//       double stress2d[3];
//       copy3Dto2Dtensors_FEEPreduced (stress[lc], stress2d);
//       errol;
//       //eqstress[lc] = sqrt( 0.5 * (stress2d));
//     }
//     else errol;
// 
//     //
//     eqstress[lc] *= this->give_volume();
//   }
// 
//   return eqstress[lc];
// }


} // end of namespace mumech

/*end of file*/