AdlibELCCoupledSolver.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) 2003-2007 California Institute of Technology
// Ralf Deiterding, ralf@amroc.net

#ifndef ADLIB_ELC_COUPLED_SOLVER_H
#define ADLIB_ELC_COUPLED_SOLVER_H

template <int dim> class AdlibCoupledSolver;

#include "external_coupling/AdlibCoupledSolver3.h"
#include "elc/LagrangianComm.h"

template <int dim>
class AdlibELCCoupledSolver : public AdlibCoupledSolver<dim> {
  typedef AdlibCoupledSolver<dim> base;
public:
  typedef elc::LagrangianComm<dim,double> elc_lagcomm_type;

  AdlibELCCoupledSolver(MPI_Comm solidComm, int numFluidNodes, int firstFluidNode) : base(solidComm) {

#ifdef DEBUG_PRINT_ELC
    ( base::log() << "*** LagrangianComm: " << numFluidNodes << "  " << firstFluidNode ).flush();
#endif  

    // instantiate an elc object for data exchange
    _elcLag = new elc_lagcomm_type(MPI_COMM_WORLD, solidComm, numFluidNodes,
                                   firstFluidNode, elc::GlobalIdentifiers);

#ifdef DEBUG_PRINT_ELC
    ( base::log() << " created.\n").flush();
#endif  
  }

  ~AdlibELCCoupledSolver() {    
    if (_elcLag) delete _elcLag;
  }
  
  virtual void sendBoundaryReceivePressure() {
    ADLIB_START_WATCH
    // update fluid/solid boundary
    adlib::AdlibBoundaryUpdate(&base::Boundary());
    
    int elNumNodes = base::Boundary().bNodes;
    int elNumElements = base::Boundary().bFacets;
    
    // list length = elNumNodes
    double *elPressures = base::Boundary().bPressure;
    int *elGlobalNodeIDs= new int[elNumNodes];
    for (int nod=0; nod<elNumNodes; ++nod) {
      int nodeNbInPart = base::Boundary().bNodeMap[nod];
      elGlobalNodeIDs[nod] =adlib::OTArray[nodeNbInPart-1]->GlobalId;
    }

    // list length = dim*elNumNodes 
    double *elCoordinates = new double[elNumNodes*dim]; ;
    computeBoundaryCoordinates(base::Boundary().bCoordinates, elCoordinates, 
                               base::Boundary().bNodeMap, elNumNodes);

    double *elVelocities = base::Boundary().bVelocities;
    
    // list length = dim*elNumElements
    int *elConnectivity = new int[dim*elNumElements]; 
    int *bIds           = base::Boundary().bConnectivities;
    int *bNodeMap       = base::Boundary().bNodeMap; 
    int *gIds           = elConnectivity; 

    for (int elem=0; elem<elNumElements; ++elem, bIds+=dim, gIds+=dim) 
      for (int a=0; a<dim; ++a) {
        int nodeNbInPart = bNodeMap[bIds[a]-1];
        gIds[a] = adlib::OTArray[nodeNbInPart-1]->GlobalId;
      }
    ADLIB_END_WATCH(CPL_SEND_OVERHEAD)
         
#ifdef DEBUG_PRINT_ELC
    ( base::log() << "*** AdlibCoupledSolver::sendBoundaryReceivePressure" ).flush();
#endif
        
    ModifyELCSendData(_elcLag);

    ADLIB_START_WATCH
      _elcLag->sendMesh(elNumNodes, reinterpret_cast<void*>(elGlobalNodeIDs), 
                        reinterpret_cast<void*>(elCoordinates), 
                        reinterpret_cast<void*>(elVelocities), 
                        elNumElements, reinterpret_cast<void*>(elConnectivity));
      _elcLag->waitForMesh();
      _elcLag->receivePressure(elNumNodes, reinterpret_cast<void*>(elPressures));
      _elcLag->waitForPressure();
    ADLIB_END_WATCH(ELC_SNDBND_RCVPRS)

    ModifyELCReceiveData(_elcLag);

    ADLIB_START_WATCH
    // Necessary if solid domain larger than fluid domain
    for (register int n=0; n<elNumNodes; n++) 
      if (elPressures[n] == std::numeric_limits<double>::max()) 
        elPressures[n] = 0.0;
    
    std::memset(adlib::forces, 0,
                sizeof(double)*(adlib::nodes)*(adlib::dof_node));

    // apply pressure loading
    adlib::AdlibBoundaryApplyPressure(&base::Boundary(),adlib::solidFluid);
    
    delete [] elConnectivity;
    delete [] elGlobalNodeIDs;
    delete [] elCoordinates;

    ADLIB_END_WATCH(CPL_RECEIVE_OVERHEAD)

#ifdef DEBUG_PRINT_ELC
    ( base::log() << " done.\n" ).flush();
#endif
  }

  virtual void computeBoundaryCoordinates(double *original, double *final, 
                                          int *map, int size)
  {
    for(int i = 0; i < size; i++)
      {
        double x0 = adlib::coordinates[adlib::dof_node*(map[i]-1)];
        double y0 = adlib::coordinates[adlib::dof_node*(map[i]-1)+1];
        double z0 = adlib::coordinates[adlib::dof_node*(map[i]-1)+2];
        double x = original[dim*i+0];
        double y = original[dim*i+1];
        double z = original[dim*i+2];
        
        double dx = x0 - base::xCenter;
        if(dx>0.)
          dx *= base::ratioXPos;
        else
          dx *= base::ratioXNeg;
        final[dim*i+0] = base::xCenter + dx + (x-x0);

        double dy = y0 - base::yCenter;
        if(dy>0.)
          dy *= base::ratioYPos;
        else
          dy *= base::ratioYNeg;
        final[dim*i+1] = base::yCenter + dy + (y-y0);

        double dz = z0 - base::zCenter;
        if(dz>0.)
          dz *= base::ratioZPos;
        else
          dz *= base::ratioZNeg;
        final[dim*i+2] = base::zCenter + dz + (z-z0);
      }
  }

  virtual void ModifyELCReceiveData(elc_lagcomm_type* elcLag) {}
  virtual void ModifyELCSendData(elc_lagcomm_type* elcLag) {}

protected:
  elc_lagcomm_type* _elcLag;
};

#endif

---