src/base/stressbalance/FETools.cc

Go to the documentation of this file.

// Copyright (C) 2009--2011 Jed Brown and Ed Bueler and Constantine Khroulev and David Maxwell
//
// This file is part of PISM.
//
// PISM is free software; you can redistribute it and/or modify it under the
// terms of the GNU General Public License as published by the Free Software
// Foundation; either version 2 of the License, or (at your option) any later
// version.
//
// PISM 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 General Public License for more
// details.
//
// You should have received a copy of the GNU General Public License
// along with PISM; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

// Utility functions used by the SSAFEM code.

#include "FETools.hh"

const FEShapeQ1::ShapeFunctionSpec FEShapeQ1::shapeFunction[FEShapeQ1::Nk] = 
{FEShapeQ1::shape0, FEShapeQ1::shape1, FEShapeQ1::shape2, FEShapeQ1::shape3};


FEElementMap::FEElementMap(const IceGrid &g)
{
  // Start by assuming ghost elements exist in all directions.
  // Elements are indexed by their lower left vertex.  If there is a ghost
  // element on the right, its i-index will be the same as the maximum
  // i-index of a non-ghost vertex in the local grid.
  xs= g.xs-1;                    // Start at ghost to the left. 
  PetscInt xf = g.xs + g.xm - 1; // End at ghost to the right.
  ys= g.ys-1;                    // Start at ghost at the bottom.
  PetscInt yf = g.ys + g.ym - 1; // End at ghost at the top.

  lxs = g.xs;
  PetscInt lxf = lxs + g.xm - 1;
  lys = g.ys;
  PetscInt lyf = lys + g.ym - 1;

  // Now correct if needed. The only way there will not be ghosts is if the 
  // grid is not periodic and we are up against the grid boundary.
  
  if( !(g.periodicity & X_PERIODIC) )
  {
    // Leftmost element has x-index 0.
    if(xs < 0){
      xs = 0;
    }
    // Rightmost vertex has index g.Mx-1, so the rightmost element has index g.Mx-2
    if(xf > g.Mx-2)
    {
      xf  = g.Mx-2;
      lxf = g.Mx-2;
    }
  }

  if( !(g.periodicity & Y_PERIODIC) )
  {
    // Bottom element has y-index 0.
    if(ys < 0){
      ys = 0;
    }
    // Topmost vertex has index g.My-1, so the topmost element has index g.My-2
    if(yf > g.My-2)
    {
      yf  = g.My-2;
      lyf = g.My-2;
    }
  }
  
  // Tally up the number of elements in each direction
  xm = xf-xs+1;
  ym = yf-ys+1;
  lxm = lxf-lxs+1;
  lym = lyf-lys+1;

}


void FEDOFMap::extractLocalDOFs(PetscInt i,PetscInt j, PetscReal const*const*xg,PetscReal *x) const
{
  x[0] = xg[i][j]; x[1] = xg[i+1][j]; x[2] = xg[i+1][j+1]; x[3] = xg[i][j+1];  
}

void FEDOFMap::extractLocalDOFs(PetscInt i,PetscInt j, PISMVector2 const*const*xg,PISMVector2 *x) const
{
  x[0] = xg[i][j]; x[1] = xg[i+1][j]; x[2] = xg[i+1][j+1]; x[3] = xg[i][j+1];    
}


void FEDOFMap::extractLocalDOFs(PetscReal const*const*xg,PetscReal *x) const
{
  extractLocalDOFs(m_i,m_j,xg,x);
}

void FEDOFMap::extractLocalDOFs(PISMVector2 const*const*xg,PISMVector2 *x) const
{
  extractLocalDOFs(m_i,m_j,xg,x);
}

void FEDOFMap::localToGlobal(PetscInt k, PetscInt *i, PetscInt *j)
{
  *i = m_i + kIOffset[k];
  *j = m_j + kJOffset[k];  
}

void FEDOFMap::reset(PetscInt i, PetscInt j, const IceGrid &grid)
{
  m_i = i; m_j = j;
  // The meaning of i and j for a PISM IceGrid and for a Petsc DA are swapped (the so-called
  // fundamental transpose.  The interface between PISM and Petsc is the stencils, so all
  // interactions with the stencils involve a transpose.
  m_col[0].j = i;   m_col[0].i = j;
  m_col[1].j = i+1; m_col[1].i = j;
  m_col[2].j = i+1; m_col[2].i = j+1;
  m_col[3].j = i;   m_col[3].i = j+1;

  memcpy(m_row,m_col,Nk*sizeof(m_col[0]));

  // We do not ever sum into rows that are not owned by the local rank.
  for(PetscInt k=0; k<Nk; k++) {         
    PetscInt pism_i = m_row[k].j, pism_j = m_row[k].i;
    if (  pism_i < grid.xs || grid.xs+grid.xm-1 < pism_i || 
                pism_j < grid.ys || grid.ys+grid.ym-1 < pism_j ) {
      markRowInvalid(k);      
    }
  }
}

void FEDOFMap::markRowInvalid(PetscInt k)
{
  m_row[k].i=m_row[k].j = kDofInvalid;
}

void FEDOFMap::markColInvalid(PetscInt k)
{
  m_col[k].i=m_col[k].j = kDofInvalid;
}

void FEDOFMap::addLocalResidualBlock(const PISMVector2 *y, PISMVector2 **yg)
{
  for (int k=0; k<Nk; k++) {
    if (m_row[k].i == kDofInvalid || m_row[k].j == kDofInvalid) continue;
    yg[m_row[k].j][m_row[k].i].u += y[k].u;
    yg[m_row[k].j][m_row[k].i].v += y[k].v;
  }
}
void FEDOFMap::addLocalResidualBlock(const PetscScalar *y, PetscScalar **yg)
{
  for (int k=0; k<Nk; k++) {
    if (m_row[k].i == kDofInvalid || m_row[k].j == kDofInvalid) continue;
    yg[m_row[k].j][m_row[k].i] += y[k];
  }
}


PetscErrorCode FEDOFMap::addLocalJacobianBlock(const PetscReal *K, Mat J)
{
  PetscErrorCode ierr = MatSetValuesBlockedStencil(J,Nk,m_row,Nk,m_col,K,ADD_VALUES);CHKERRQ(ierr);  
  return 0;
}


PetscErrorCode FEDOFMap::setJacobianDiag(PetscInt i, PetscInt j, const PetscReal*K, Mat J)
{
  MatStencil row;
  row.i=j; row.j=i;
  PetscErrorCode ierr = MatSetValuesBlockedStencil(J,1,&row,1,&row,K,INSERT_VALUES);CHKERRQ(ierr);  
  return 0;
}

const PetscInt FEDOFMap::kIOffset[4] = {0,1,1,0};
const PetscInt FEDOFMap::kJOffset[4] = {0,0,1,1};


FEQuadrature::FEQuadrature()
{  
}

void FEQuadrature::getWeightedJacobian(PetscReal *jxw)
{
  for(int q=0;q<Nq;q++)
  {
    jxw[q] = m_jacobianDet * quadWeights[q];
  } 
}

void FEQuadrature::init(const IceGrid &grid,PetscScalar L)
{
  // Since we use uniform cartesian coordinates, the Jacobian is constant and diagonal on every element.
  // Note that the reference element is \f$ [-1,1]^2 \f$ hence the extra factor of 1/2.
  PetscReal jacobian_x = 0.5*grid.dx/L;
  PetscReal jacobian_y = 0.5*grid.dy/L;
  m_jacobianDet = jacobian_x*jacobian_y;

  FEShapeQ1 shape;
  for(int q=0; q<Nq; q++){
    for(int k=0; k<Nk; k++){
      shape.eval(k,quadPoints[q][0],quadPoints[q][1],&m_germs[q][k]);
      m_germs[q][k].dx /= jacobian_x;
      m_germs[q][k].dy /= jacobian_y;
    }
  }
}

//* The return value is an Nq by Nk array of FEFunctionGerms. */
const FEFunctionGerm (*FEQuadrature::testFunctionValues())[FEQuadrature::Nq]
{
  return m_germs;
}

//* The return value is an array of Nk FEFunctionGerms. */
const FEFunctionGerm *FEQuadrature::testFunctionValues(PetscInt q)
{
  return m_germs[q];
}

const FEFunctionGerm *FEQuadrature::testFunctionValues(PetscInt q, PetscInt k)
{
  return m_germs[q] + k;
}


void FEQuadrature::computeTrialFunctionValues(const PetscReal *x, PetscReal *vals)
{
  for (int q=0; q<Nq; q++) {
    const FEFunctionGerm *test = m_germs[q];
    vals[q] = 0;
    for (int k=0; k<Nk; k++) {
      vals[q] += test[k].val * x[k];
    }
  }
}

void FEQuadrature::computeTrialFunctionValues(const PetscReal *x, PetscReal *vals, PetscReal *dx, PetscReal *dy)
{
  for (int q=0; q<Nq; q++) {
    const FEFunctionGerm *test = m_germs[q];
    vals[q] = 0; dx[q] = 0; dy[q] = 0;
    for (int k=0; k<Nk; k++) {
      vals[q] += test[k].val * x[k];
      dx[q]   += test[k].dx * x[k];
      dy[q]   += test[k].dy * x[k];
    }
  }
}

void FEQuadrature::computeTrialFunctionValues( PetscInt i, PetscInt j, const FEDOFMap &dof, 
                                             PetscReal const*const*xg, PetscReal *vals)
{
  dof.extractLocalDOFs(i,j,xg,m_tmpScalar);
  computeTrialFunctionValues(m_tmpScalar,vals);
}

void FEQuadrature::computeTrialFunctionValues( PetscInt i, PetscInt j, const FEDOFMap &dof, PetscReal const*const*xg, 
                                 PetscReal *vals, PetscReal *dx, PetscReal *dy)
{
  dof.extractLocalDOFs(i,j,xg,m_tmpScalar);
  computeTrialFunctionValues(m_tmpScalar,vals,dx,dy);  
}

void FEQuadrature::computeTrialFunctionValues( const PISMVector2 *x,  PISMVector2 *vals)
{
  for (int q=0; q<Nq; q++) {
    vals[q].u = 0; vals[q].v = 0;
    const FEFunctionGerm *test = m_germs[q];
    for (int k=0; k<Nk; k++) {
      vals[q].u += test[k].val * x[k].u;
      vals[q].v += test[k].val * x[k].v;
    }
  }  

}

void FEQuadrature::computeTrialFunctionValues( const PISMVector2 *x,  PISMVector2 *vals, PetscReal (*Dv)[3] )
{
  for (int q=0; q<Nq; q++) {
    vals[q].u = 0; vals[q].v = 0;
    PetscReal *Dvq = Dv[q];
    Dvq[0]=0; Dvq[1]=0; Dvq[2]=0;
    const FEFunctionGerm *test = m_germs[q];
    for(int k=0; k<Nk; k++) {
      vals[q].u += test[k].val * x[k].u;
      vals[q].v += test[k].val * x[k].v;
      Dvq[0] += test[k].dx * x[k].u;
      Dvq[1] += test[k].dy * x[k].v;
      Dvq[2] += 0.5*(test[k].dy*x[k].u + test[k].dx*x[k].v);
    }
  }  
}

void FEQuadrature::computeTrialFunctionValues( const PISMVector2 *x,  PISMVector2 *vals, PISMVector2 *dx, PISMVector2 *dy)
{
  for (int q=0; q<Nq; q++) {
    vals[q].u = 0; vals[q].v = 0;
    dx[q].u = 0; dx[q].v = 0;
    dy[q].u = 0; dy[q].v = 0;
    const FEFunctionGerm *test = m_germs[q];
    for(int k=0; k<Nk; k++) {
      vals[q].u += test[k].val * x[k].u;
      vals[q].v += test[k].val * x[k].v;
      dx[q].u += test[k].dx * x[k].u;
      dx[q].v += test[k].dx * x[k].v;
      dy[q].u += test[k].dy * x[k].u;
      dy[q].v += test[k].dy * x[k].v;
    }
  }    
}
  

void FEQuadrature::computeTrialFunctionValues( PetscInt i, PetscInt j, const FEDOFMap &dof,                                              
                                             PISMVector2 const*const*xg, PISMVector2 *vals )
{
  dof.extractLocalDOFs(i,j,xg,m_tmpVector);
  computeTrialFunctionValues(m_tmpVector,vals);
}

void FEQuadrature::computeTrialFunctionValues( PetscInt i, PetscInt j, const FEDOFMap &dof,                                              
                                             PISMVector2 const*const*xg, PISMVector2 *vals, PetscReal (*Dv)[3] )
{
  dof.extractLocalDOFs(i,j,xg,m_tmpVector);
  computeTrialFunctionValues(m_tmpVector,vals,Dv);
}

const PetscReal FEQuadrature::quadPoints[FEQuadrature::Nq][2] = 
                                          {{ -0.57735026918962573, -0.57735026918962573 },
                                           {  0.57735026918962573, -0.57735026918962573 },
                                           {  0.57735026918962573,  0.57735026918962573 },
                                           { -0.57735026918962573,  0.57735026918962573 }};

const PetscReal FEQuadrature::quadWeights[FEQuadrature::Nq]  = {1,1,1,1};



PetscTruth Floating(const IceFlowLaw &ice, PetscScalar ocean_rho,
                           PetscReal H, PetscReal bed)
{
  return ice.rho*H + ocean_rho*bed < 0 ? PETSC_TRUE : PETSC_FALSE;
}


int PismIntMask(PetscScalar maskvalue) {
  return static_cast<int>(floor(maskvalue + 0.5));
}