<commit_before>/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    itkFEMElementMembraneC02D.cxx
  Language:  C++
  Date:      $Date$
  Version:   $Revision$

  Copyright (c) 2002 Insight Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/

// disable debug warnings in MS compiler
#ifdef _MSC_VER
#pragma warning(disable: 4786)
#endif

#include "itkFEMElementMembraneC02D.h"
#include "FEM/itkFEMLoadGrav.h"
#include "FEM/itkFEMLoadEdge.h"
#include "FEM/itkFEMObjectFactory.h"
#include "FEM/itkFEMUtility.h"
#include "vnl/vnl_math.h"

namespace itk {
namespace fem {




/**
 * Constructor for class MembraneC02D
 */
MembraneC02D::MembraneC02D(  Node::ConstPointer n1_,
          Node::ConstPointer n2_,
          Node::ConstPointer n3_,
          Node::ConstPointer n4_,
          Material::ConstPointer p_)
{
  try
  {
    /**
     * Initialize the pointers to nodes and check that
     * we were given the pointers to the right node class.
     * if the node class was incorrect a bad_cast exception is thrown
     */
    m_node1=&dynamic_cast<const NodeXY&>(*n1_);
    m_node2=&dynamic_cast<const NodeXY&>(*n2_);
    m_node3=&dynamic_cast<const NodeXY&>(*n3_);
    m_node4=&dynamic_cast<const NodeXY&>(*n4_);
    m_mat=&dynamic_cast<const MaterialStandard&>(*p_);
  }
  catch ( std::bad_cast )
  {
    throw FEMExceptionWrongClass(__FILE__,__LINE__,"MembraneC02D::MembraneC02D()");
  }

}




/**
 * Returns the stiffness matrix for MembraneC02D element
 */
vnl_matrix<MembraneC02D::Float> MembraneC02D::Ke() const 
{
  vnl_matrix<Float> MatKe(8,8), I(2,2), shapeD(4,2), shapeINVD(4,2),
            J(2,2), D(3,3,0.0), B(3,8), DB(3,8);

  Float detJ;
  
  /** Gaussian integration points */
  Float pt = 1.0 / sqrt(3.0);
  Float GPoints[][2] = {{-pt, -pt}, {pt, -pt}, {pt, pt}, {-pt, pt}};

  /** Material properties matrix */
  Float disot = m_mat->E;
    
  D(0,0) = disot;
  D(1,1) = disot;
  D(2,2) = disot;

  Float x[2];

  /** Initialize stiffness matrix */
  MatKe.fill(0.0);

  /** For each integration point */
  for (int k=0; k<4; k++) {
    /** Get the k-th integration point */
    x[0] = GPoints[k][0];
    x[1] = GPoints[k][1];


    /**
     * Computes the Jacobian matrix and its determinant
     * at the k-th integration point
     */
    J = ComputeJacobianMatrixAt(x);
    detJ = JacobianMatrixDeterminant(J);

    /** Computes the inverse of the Jacobian matrix */
    I = ComputeJacobianInverse(J, detJ);

    /** Computes the shape function derivatives at integration point */
    shapeD = ComputeShapeFunctionDerivativesAt(x);
    /**
     * Computes the shape function derivatives in Cartesian coordinates
     * at integration point
     */
    shapeINVD = ComputeShapeFunctionCartDerivatives(I, shapeD);

    /** Computes the strain (B) matrix */
    B = ComputeBMatrix(shapeINVD);
//  std::cout << B << std::endl;
    /** Computes the matrix multiplication DB */
 //   DB = ComputeDBMatrix(D,B);

    /**
     * Add the contribution of k-th integration point to
     * the stiffness matrix
     */
  MatKe += detJ*(B.transpose()*(D*B));
  }
  return MatKe;
}



/**
 * Returns the stiffness matrix for MembraneC02D element
 */
vnl_matrix<MembraneC02D::Float> MembraneC02D::Me() const 
{
  vnl_matrix<Float> MatMe(8,8), I(2,2),  shapeFINV(4,2), 
     shapeF(4,2), J(2,2), Nmat(3,8);

  Float detJ,rho=1.;
  
  /** Gaussian integration points */
  Float pt = 1.0 / sqrt(3.0);
  Float GPoints[][2] = {{-pt, -pt}, {pt, -pt}, {pt, pt}, {-pt, pt}};


  Float x[2];

  /** Initialize stiffness matrix */
  MatMe.fill(0.0);

  /** For each integration point */
  for (int k=0; k<4; k++) {
    /** Get the k-th integration point */
    x[0] = GPoints[k][0];
    x[1] = GPoints[k][1];


    /**
     * Computes the Jacobian matrix and its determinant
     * at the k-th integration point
     */
    J = ComputeJacobianMatrixAt(x);
    detJ = JacobianMatrixDeterminant(J);
<commit_msg>Commented unused variable rho<commit_after>/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    itkFEMElementMembraneC02D.cxx
  Language:  C++
  Date:      $Date$
  Version:   $Revision$

  Copyright (c) 2002 Insight Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/

// disable debug warnings in MS compiler
#ifdef _MSC_VER
#pragma warning(disable: 4786)
#endif

#include "itkFEMElementMembraneC02D.h"
#include "FEM/itkFEMLoadGrav.h"
#include "FEM/itkFEMLoadEdge.h"
#include "FEM/itkFEMObjectFactory.h"
#include "FEM/itkFEMUtility.h"
#include "vnl/vnl_math.h"

namespace itk {
namespace fem {




/**
 * Constructor for class MembraneC02D
 */
MembraneC02D::MembraneC02D(  Node::ConstPointer n1_,
          Node::ConstPointer n2_,
          Node::ConstPointer n3_,
          Node::ConstPointer n4_,
          Material::ConstPointer p_)
{
  try
  {
    /**
     * Initialize the pointers to nodes and check that
     * we were given the pointers to the right node class.
     * if the node class was incorrect a bad_cast exception is thrown
     */
    m_node1=&dynamic_cast<const NodeXY&>(*n1_);
    m_node2=&dynamic_cast<const NodeXY&>(*n2_);
    m_node3=&dynamic_cast<const NodeXY&>(*n3_);
    m_node4=&dynamic_cast<const NodeXY&>(*n4_);
    m_mat=&dynamic_cast<const MaterialStandard&>(*p_);
  }
  catch ( std::bad_cast )
  {
    throw FEMExceptionWrongClass(__FILE__,__LINE__,"MembraneC02D::MembraneC02D()");
  }

}




/**
 * Returns the stiffness matrix for MembraneC02D element
 */
vnl_matrix<MembraneC02D::Float> MembraneC02D::Ke() const 
{
  vnl_matrix<Float> MatKe(8,8), I(2,2), shapeD(4,2), shapeINVD(4,2),
            J(2,2), D(3,3,0.0), B(3,8), DB(3,8);

  Float detJ;
  
  /** Gaussian integration points */
  Float pt = 1.0 / sqrt(3.0);
  Float GPoints[][2] = {{-pt, -pt}, {pt, -pt}, {pt, pt}, {-pt, pt}};

  /** Material properties matrix */
  Float disot = m_mat->E;
    
  D(0,0) = disot;
  D(1,1) = disot;
  D(2,2) = disot;

  Float x[2];

  /** Initialize stiffness matrix */
  MatKe.fill(0.0);

  /** For each integration point */
  for (int k=0; k<4; k++) {
    /** Get the k-th integration point */
    x[0] = GPoints[k][0];
    x[1] = GPoints[k][1];


    /**
     * Computes the Jacobian matrix and its determinant
     * at the k-th integration point
     */
    J = ComputeJacobianMatrixAt(x);
    detJ = JacobianMatrixDeterminant(J);

    /** Computes the inverse of the Jacobian matrix */
    I = ComputeJacobianInverse(J, detJ);

    /** Computes the shape function derivatives at integration point */
    shapeD = ComputeShapeFunctionDerivativesAt(x);
    /**
     * Computes the shape function derivatives in Cartesian coordinates
     * at integration point
     */
    shapeINVD = ComputeShapeFunctionCartDerivatives(I, shapeD);

    /** Computes the strain (B) matrix */
    B = ComputeBMatrix(shapeINVD);
//  std::cout << B << std::endl;
    /** Computes the matrix multiplication DB */
 //   DB = ComputeDBMatrix(D,B);

    /**
     * Add the contribution of k-th integration point to
     * the stiffness matrix
     */
  MatKe += detJ*(B.transpose()*(D*B));
  }
  return MatKe;
}



/**
 * Returns the stiffness matrix for MembraneC02D element
 */
vnl_matrix<MembraneC02D::Float> MembraneC02D::Me() const 
{
  vnl_matrix<Float> MatMe(8,8), I(2,2),  shapeFINV(4,2), 
     shapeF(4,2), J(2,2), Nmat(3,8);

  Float detJ;
//  Float rho=1.;
  
  /** Gaussian integration points */
  Float pt = 1.0 / sqrt(3.0);
  Float GPoints[][2] = {{-pt, -pt}, {pt, -pt}, {pt, pt}, {-pt, pt}};


  Float x[2];

  /** Initialize stiffness matrix */
  MatMe.fill(0.0);

  /** For each integration point */
  for (int k=0; k<4; k++) {
    /** Get the k-th integration point */
    x[0] = GPoints[k][0];
    x[1] = GPoints[k][1];


    /**
     * Computes the Jacobian matrix and its determinant
     * at the k-th integration point
     */
    J = ComputeJacobianMatrixAt(x);
    detJ = JacobianMatrixDeterminant(J);

{'fraction_non_alphanumeric': 0.16951093803688008, 'fraction_numerical': 0.02542664070553201, 'mean_word_length': 3.3014778325123153, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5892357', 'n_tokens_mistral': 3159, 'n_tokens_neox': 2934, 'n_words': 755}