LAIHelperFunctions.hpp

Go to the documentation of this file.

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2024 TotalEnergies
 * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2023-2024 Chevron
 * Copyright (c) 2019-     GEOS/GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */
 
#ifndef GEOS_LINEARALGEBRA_UTILITIES_LAIHELPERFUNCTIONS_HPP_
#define GEOS_LINEARALGEBRA_UTILITIES_LAIHELPERFUNCTIONS_HPP_
 
#include "common/DataTypes.hpp"
#include "linearAlgebra/interfaces/InterfaceTypes.hpp"
#include "linearAlgebra/DofManager.hpp"
#include "mesh/MeshBody.hpp"
#include "mesh/NodeManager.hpp"
#include "mesh/ElementRegionManager.hpp"
 
namespace geos
{
namespace LAIHelperFunctions
{
 
template< typename MATRIX >
void makeIdentity( localIndex const n,
                   MPI_Comm const & comm,
                   MATRIX & mat )
{
  mat.createWithLocalSize( n, 1, comm );
  mat.open();
  for( globalIndex i = mat.ilower(); i < mat.iupper(); ++i )
  {
    mat.insert( i, i, 1.0 );
  }
  mat.close();
}
 
template< typename MATRIX >
void createPermutationMatrix( NodeManager const & nodeManager,
                              int const nDofPerNode,
                              string const & dofKey,
                              MATRIX & permutationMatrix )
{
  /* Crates a permutation matrix for a given nodal variable specified by the DofKey.
   * It consider that nDofPerNode dofs are associated with each node (e.g., nDofPerNode = 3 for the displacement).
   *
   * The permutation matrix maps from the dofs ordering set by the DOFManager to the ordering based on the global
   * indexes of the nodes.
   */
 
  localIndex const numLocalRows = nodeManager.getNumberOfLocalIndices() * nDofPerNode;
  permutationMatrix.createWithLocalSize( numLocalRows, numLocalRows, 1, MPI_COMM_GEOS );
 
  arrayView1d< globalIndex const > const & dofNumber = nodeManager.getReference< globalIndex_array >( dofKey );
  arrayView1d< globalIndex const > const & localToGlobal = nodeManager.localToGlobalMap();
 
  permutationMatrix.open();
  for( localIndex a = 0; a < nodeManager.size(); ++a )
  {
    if( dofNumber[a] >= 0 )
    {
      for( int d = 0; d < nDofPerNode; ++d )
      {
        globalIndex const rowIndex    = localToGlobal[a] * nDofPerNode + d;
        globalIndex const columnIndex = dofNumber[a] + d;
        permutationMatrix.insert( rowIndex, columnIndex, 1.0 );
      }
    }
  }
  permutationMatrix.close();
  permutationMatrix.set( 1.0 );
}
 
template< typename MATRIX >
void createPermutationMatrix( ElementRegionManager const & elemManager,
                              int const nDofPerCell,
                              string const & dofKey,
                              MATRIX & permutationMatrix )
{
  /* Crates a permutation matrix for a given cell centered variable specified by the DofKey.
   * It consider that nDofPerCell dofs are associated with each cell (e.g., nDofPerCell = 3 for the displacement).
   *
   * The permutation matrix maps from the dofs ordering set by the DOFManager to the ordering based on the global
   * indexes of the cells.
   */
 
  // Create permutation matrix
  localIndex numLocalRows = 0;
  elemManager.forElementSubRegions< ElementSubRegionBase >( [&]( ElementSubRegionBase const & elementSubRegion )
  {
    if( elementSubRegion.hasWrapper( dofKey ) )
    {
      numLocalRows += elementSubRegion.getNumberOfLocalIndices() * nDofPerCell;
    }
  } );
  permutationMatrix.createWithLocalSize( numLocalRows, numLocalRows, 1, MPI_COMM_GEOS );
 
  permutationMatrix.open();
  elemManager.forElementSubRegions< ElementSubRegionBase >( [&]( ElementSubRegionBase const & elementSubRegion )
  {
    localIndex const numElems = elementSubRegion.size();
    arrayView1d< globalIndex const > const & dofNumber = elementSubRegion.getReference< array1d< globalIndex > >( dofKey );
    arrayView1d< globalIndex const > const & localToGlobal = elementSubRegion.localToGlobalMap();
 
    for( localIndex k = 0; k < numElems; ++k )
    {
      if( dofNumber[k] >= 0 )
      {
        for( int d = 0; d < nDofPerCell; ++d )
        {
          globalIndex const rowIndex    = localToGlobal[k] * nDofPerCell + d;
          globalIndex const columnIndex = dofNumber[k] + d;
          permutationMatrix.insert( rowIndex, columnIndex, 1.0 );
        }
      }
    }
  } );
  permutationMatrix.close();
  permutationMatrix.set( 1.0 );
}
 
template< typename VECTOR, typename MATRIX >
VECTOR permuteVector( VECTOR const & vector,
                      MATRIX const & permutationMatrix )
{
  VECTOR permutedVector;
  permutedVector.create( vector.localSize(), permutationMatrix.comm() );
  permutationMatrix.apply( vector, permutedVector );
  return permutedVector;
}
 
template< typename MATRIX >
MATRIX permuteMatrix( MATRIX const & matrix,
                      MATRIX const & permutationMatrix )
{
  MATRIX permutedMatrix;
  matrix.multiplyRARt( permutationMatrix, permutedMatrix );
  return matrix;
}
 
template< typename MATRIX >
MATRIX permuteMatrix( MATRIX const & matrix,
                      MATRIX const & permutationMatrixLeft,
                      MATRIX const & permutationMatrixRight )
{
  MATRIX permutedMatrix;
  matrix.multiplyRAP( permutationMatrixLeft, permutationMatrixRight, permutedMatrix );
  return matrix;
}
 
template< typename VECTOR >
array1d< VECTOR >
computeRigidBodyModes( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition,
                       arrayView1d< globalIndex const > const & dofIndex,
                       globalIndex const dofOffset,
                       localIndex const numLocalDof )
{
  GEOS_ASSERT_EQ( nodePosition.size( 0 ), dofIndex.size() );
  integer const numComponents = nodePosition.size( 1 );
  integer const numRidigBodyModes = numComponents * ( numComponents + 1 ) / 2;
 
  array1d< VECTOR > rigidBodyModes( numRidigBodyModes );
 
  // Translation RBMs
  for( localIndex k = 0; k < numComponents; ++k )
  {
    rigidBodyModes[k].create( numLocalDof, MPI_COMM_GEOS );
    arrayView1d< real64 > const values = rigidBodyModes[k].open();
    forAll< parallelHostPolicy >( dofIndex.size(), [=]( localIndex const i )
    {
      localIndex const localDof = LvArray::integerConversion< localIndex >( dofIndex[i] - dofOffset );
      if( 0 <= localDof && localDof < numLocalDof )
      {
        values[localDof + k] = 1.0;
      }
    } );
    rigidBodyModes[k].close();
    rigidBodyModes[k].scale( 1.0 / rigidBodyModes[k].norm2() );
  }
 
  // Rotation RBMs
  for( localIndex k = numComponents; k < numRidigBodyModes; ++k )
  {
    rigidBodyModes[k].create( numLocalDof, MPI_COMM_GEOS );
    arrayView1d< real64 > const values = rigidBodyModes[k].open();
    integer const ind[2] = { ( k - numComponents + 1 ) % numComponents,
                             ( k - numComponents + 2 ) % numComponents };
    forAll< parallelHostPolicy >( dofIndex.size(), [=]( localIndex const i )
    {
      localIndex const localDof = LvArray::integerConversion< localIndex >( dofIndex[i] - dofOffset );
      if( 0 <= localDof && localDof < numLocalDof )
      {
        values[localDof + ind[0]] = -nodePosition( i, ind[1] );
        values[localDof + ind[1]] = +nodePosition( i, ind[0] );
      }
    } );
    rigidBodyModes[k].close();
    for( localIndex j = 0; j < k; ++j )
    {
      rigidBodyModes[k].axpy( -rigidBodyModes[k].dot( rigidBodyModes[j] ), rigidBodyModes[j] );
    }
    rigidBodyModes[k].scale( 1.0 / rigidBodyModes[k].norm2() );
  }
 
  return rigidBodyModes;
}
 
} // LAIHelperFunctions namespace
 
} // geos namespace
 
#endif /*GEOS_LINEARALGEBRA_UTILITIES_LAIHELPERFUNCTIONS_HPP_*/