isotropicelasticeqvar.cpp

Go to the documentation of this file.

#include "isotropicelasticeqvar.h"


const Index IsotropicElasticEqVar::m_local_dof_face_cell_map[6][4] = {{0,2,4,6},
                                                                       {1,3,5,7},
                                                                       {0,1,4,5},
                                                                       {2,3,6,7},
                                                                       {0,1,2,3},
                                                                       {4,5,6,7}};




IsotropicElasticEqVar::IsotropicElasticEqVar(OrthoMesh &mesh,  VecDouble &cK, double l1,double l2, Function3D &fD, Function3D &fN,unsigned debugLevel) 
   :m_FE(1),m_FES(m_FE,3),m_cQG(2),m_fQG(2), //Initialize finite element and quadrature rules for cells and faces
    m_global_dofs_map(OrthoMesh::VERTICES_PER_CELL*3),m_mesh(mesh),m_l1(l1),m_l2(l2)
 {
   initFEValues(m_mesh,m_FES,m_cQG,m_fQG, update_quadrature_points|update_gradients|update_JxW_values|update_values, update_JxW_values|update_values|update_quadrature_points);

 }








// Index IsotropicElasticEqVar::local_face_to_local_cell_map(Index face,Index local_dof_face,unsigned component) 
// {
//   assert(component<3);
//   assert(face < OrthoMesh::FACES_PER_CELL);
//   assert(local_dof_face<n_dofs_per_face());
//   div_t _div = div(local_dof_face,4);
//   return m_local_dof_face_cell_map[face][_div.rem] + _div.quot*8;

// }


void IsotropicElasticEqVar::setCell(OrthoMesh::Cell_It &cell) 
{
  DealFEWrapper::setCell(cell);
  unsigned i=0;
  //Make the local to global dof mapping
    for (unsigned v=0;v<OrthoMesh::VERTICES_PER_CELL;v++)
    {
      for (unsigned comp=0;comp<3;comp++)
      {
        m_global_dofs_map[i++]=cell->vertex_index(static_cast<VertexDirection3D>(v))+m_mesh.numVertices()*comp;
      }
    }
}


double IsotropicElasticEqVar::localEqVar(unsigned i,unsigned j,unsigned qPoint) 
{
  unsigned iComp=local_dof_component(i);
  unsigned jComp=local_dof_component(j);
  const  Tensor<1,DIM>& gradV = shape_grad(i,qPoint);
  double l1 = m_l1;
  double l2 = m_l2;
  switch(iComp)
  {
  case U0:
    {

      switch(jComp)
      {
      case U0:
        return (
                (l1+2*l2)*shape_grad(j,qPoint)[0]*gradV[0] +
                l2*shape_grad(j,qPoint)[1]*gradV[1] +
                l2*shape_grad(j,qPoint)[2]*gradV[2])*JxW(qPoint);
      case U1:
        return (
                l1*shape_grad(j,qPoint)[1]*gradV[0]+
                l2*shape_grad(j,qPoint)[0]*gradV[1])*JxW(qPoint);
      case U2:
        return (
                l1*shape_grad(j,qPoint)[2]*gradV[0] +
                l2*shape_grad(j,qPoint)[0]*gradV[2])*JxW(qPoint);
      }
    }
  case U1:
    {
      switch (jComp)
      {
      case U0:
        return (
                l2*shape_grad(j,qPoint)[1]*gradV[0] +
                l1*shape_grad(j,qPoint)[0]*gradV[1]
                )*JxW(qPoint); 
      case U1:
        return (
                l2       *shape_grad(j,qPoint)[0]*gradV[0]+
                (l1+2*l2)*shape_grad(j,qPoint)[1]*gradV[1]+
                l2       *shape_grad(j,qPoint)[2]*gradV[2]
                )*JxW(qPoint);
      case U2:
        return (
                l1*shape_grad(j,qPoint)[2]*gradV[1] +
                l2*shape_grad(j,qPoint)[1]*gradV[2]
                )*JxW(qPoint);
      }
    }
  case U2: //To do 
    {
      switch (jComp)
      {
      case U0:
        return (
                l2*shape_grad(j,qPoint)[2]*gradV[0] +
                l1*shape_grad(j,qPoint)[0]*gradV[2]
                )*JxW(qPoint); 
      case U1:
        return ( 
                l2*shape_grad(j,qPoint)[2]*gradV[1] +
                l1*shape_grad(j,qPoint)[1]*gradV[2]
                 )*JxW(qPoint);
      case U2:
        return (
                l2       *shape_grad(j,qPoint)[0]*gradV[0]+
                l2       *shape_grad(j,qPoint)[1]*gradV[1]+
                (l1+2*l2)*shape_grad(j,qPoint)[2]*gradV[2]
                )*JxW(qPoint);
      }
    }
  }
  assert(0); //If reached this part there is a bug.
  return NAN;
}




int IsotropicElasticEqVar::nullProduct(unsigned i,unsigned j) 
{
return false;
}


double IsotropicElasticEqVar::localRHS(unsigned i,unsigned j) 
{
return 0.0;
}

---