void mitk::PolyDataGLMapper2D::Paint( mitk::BaseRenderer * renderer )
{
    if ( IsVisible( renderer ) == false )
        return ;

    // ok, das ist aus GenerateData kopiert
    mitk::BaseData::Pointer input = const_cast<mitk::BaseData*>( GetData() );

    assert( input );

    input->Update();

    vtkPolyData * vtkpolydata = this->GetVtkPolyData();
    assert( vtkpolydata );


    vtkLinearTransform * vtktransform = GetDataNode() ->GetVtkTransform();

    if (vtktransform)
    {
      vtkLinearTransform * inversetransform = vtktransform->GetLinearInverse();

      Geometry2D::ConstPointer worldGeometry = renderer->GetCurrentWorldGeometry2D();
      PlaneGeometry::ConstPointer worldPlaneGeometry = dynamic_cast<const PlaneGeometry*>( worldGeometry.GetPointer() );

      if ( vtkpolydata != NULL )
      {
          Point3D point;
          Vector3D normal;

          if(worldPlaneGeometry.IsNotNull())
          {
            // set up vtkPlane according to worldGeometry
            point=worldPlaneGeometry->GetOrigin();
            normal=worldPlaneGeometry->GetNormal(); normal.Normalize();
            m_Plane->SetTransform((vtkAbstractTransform*)NULL);
          }
          else
          {
            //@FIXME: does not work correctly. Does m_Plane->SetTransform really transforms a "plane plane" into a "curved plane"?
            return;
            AbstractTransformGeometry::ConstPointer worldAbstractGeometry = dynamic_cast<const AbstractTransformGeometry*>(renderer->GetCurrentWorldGeometry2D());
            if(worldAbstractGeometry.IsNotNull())
            {
              // set up vtkPlane according to worldGeometry
              point=const_cast<mitk::BoundingBox*>(worldAbstractGeometry->GetParametricBoundingBox())->GetMinimum();
              FillVector3D(normal, 0, 0, 1);
              m_Plane->SetTransform(worldAbstractGeometry->GetVtkAbstractTransform()->GetInverse());
            }
            else
              return;
          }

          vtkFloatingPointType vp[ 3 ], vnormal[ 3 ];

          vnl2vtk(point.Get_vnl_vector(), vp);
          vnl2vtk(normal.Get_vnl_vector(), vnormal);

          //normally, we would need to transform the surface and cut the transformed surface with the cutter.
          //This might be quite slow. Thus, the idea is, to perform an inverse transform of the plane instead.
          //@todo It probably does not work for scaling operations yet:scaling operations have to be
          //dealed with after the cut is performed by scaling the contour.
          inversetransform->TransformPoint( vp, vp );
          inversetransform->TransformNormalAtPoint( vp, vnormal, vnormal );

          m_Plane->SetOrigin( vp );
          m_Plane->SetNormal( vnormal );

          // set data into cutter
          m_Cutter->SetInput( vtkpolydata );
          //    m_Cutter->GenerateCutScalarsOff();
          //    m_Cutter->SetSortByToSortByCell();

          // calculate the cut
          m_Cutter->Update();

          // fetch geometry
          mitk::DisplayGeometry::Pointer displayGeometry = renderer->GetDisplayGeometry();
          assert( displayGeometry );
          //  float toGL=displayGeometry->GetSizeInDisplayUnits()[1];

          //apply color and opacity read from the PropertyList
          ApplyProperties( renderer );

          // traverse the cut contour
          vtkPolyData * contour = m_Cutter->GetOutput();

          vtkPoints *vpoints = contour->GetPoints();
          vtkCellArray *vpolys = contour->GetLines();
          vtkPointData *vpointdata = contour->GetPointData();
          vtkDataArray* vscalars = vpointdata->GetScalars();

          vtkCellData *vcelldata = contour->GetCellData();
          vtkDataArray* vcellscalars = vcelldata->GetScalars();

          int i, numberOfCells = vpolys->GetNumberOfCells();

          Point3D p;
          Point2D p2d, last, first;

//.........这里部分代码省略.........

template<> void AddUnitState(Creature* obj, CellCoord const& cellCoord)
{
    Cell cell(cellCoord);

    obj->SetCurrentCell(cell);
}

void mitk::SurfaceGLMapper2D::PaintCells(mitk::BaseRenderer* renderer, vtkPolyData* contour,
                                       const PlaneGeometry* worldGeometry,
                                       const DisplayGeometry* displayGeometry,
                                       vtkLinearTransform * vtktransform,
                                       vtkLookupTable *lut,
                                       vtkPolyData* original3DObject)
{
  // deprecated settings
  bool usePointData = false;

  bool useCellData = false;
  this->GetDataNode()->GetBoolProperty("deprecated useCellDataForColouring", useCellData, renderer);

  bool scalarVisibility = false;
  this->GetDataNode()->GetBoolProperty("scalar visibility", scalarVisibility, renderer);

  if(scalarVisibility)
  {
    VtkScalarModeProperty* scalarMode;
    if(this->GetDataNode()->GetProperty(scalarMode, "scalar mode", renderer))
    {
      if( (scalarMode->GetVtkScalarMode() == VTK_SCALAR_MODE_USE_POINT_DATA) ||
        (scalarMode->GetVtkScalarMode() == VTK_SCALAR_MODE_DEFAULT) )
      {
        usePointData = true;
      }
      if(scalarMode->GetVtkScalarMode() == VTK_SCALAR_MODE_USE_CELL_DATA)
      {
        useCellData = true;
      }
    }
    else
    {
      usePointData = true;
    }
  }

  vtkPoints    *vpoints = contour->GetPoints();
  vtkDataArray *vpointscalars = contour->GetPointData()->GetScalars();

  vtkCellArray *vlines  = contour->GetLines();
  vtkDataArray* vcellscalars = contour->GetCellData()->GetScalars();

  Point3D p; Point2D p2d, last;
  int i, j;
  int numberOfLines = vlines->GetNumberOfCells();

  glLineWidth( m_LineWidth );
  glBegin (GL_LINES);

  glColor4fv(m_LineColor);

  double distanceSinceLastNormal(0.0);

  vlines->InitTraversal();
  for(i=0;i<numberOfLines;++i)
  {
    vtkIdType *cell(NULL);
    vtkIdType cellSize(0);
    double vp[3];

    vlines->GetNextCell(cellSize, cell);

    vpoints->GetPoint(cell[0], vp);
    //take transformation via vtktransform into account
    vtktransform->TransformPoint(vp, vp);
    vtk2itk(vp, p);

    //convert 3D point (in mm) to 2D point on slice (also in mm)
    worldGeometry->Map(p, p2d);

    //convert point (until now mm and in world coordinates) to display coordinates (units )
    displayGeometry->WorldToDisplay(p2d, p2d);
    last=p2d;

    for(j=1; j<cellSize; ++j)
    {
      vpoints->GetPoint(cell[j], vp);
      Point3D originalPoint;
      vtk2itk(vp, originalPoint);
      //take transformation via vtktransform into account
      vtktransform->TransformPoint(vp, vp);
      vtk2itk(vp, p);

      //convert 3D point (in mm) to 2D point on slice (also in mm)
      worldGeometry->Map(p, p2d);

      //convert point (until now mm and in world coordinates) to display coordinates (units )
      displayGeometry->WorldToDisplay(p2d, p2d);

      double color[3];
      if (useCellData && vcellscalars != NULL )
      {
        // color each cell according to cell data
        lut->GetColor( vcellscalars->GetComponent(i,0),color);
        glColor3f(color[0],color[1],color[2]);
        glVertex2f(last[0], last[1]);
        glVertex2f(p2d[0], p2d[1]);
      }
      else if (usePointData && vpointscalars != NULL )
//.........这里部分代码省略.........

    static bool HandleGPSCommand(ChatHandler* handler, char const* args)
    {
        WorldObject* object = NULL;
        if (*args)
        {
            uint64 guid = handler->extractGuidFromLink((char*)args);
            if (guid)
                object = (WorldObject*)ObjectAccessor::GetObjectByTypeMask(*handler->GetSession()->GetPlayer(), guid, TYPEMASK_UNIT | TYPEMASK_GAMEOBJECT);

            if (!object)
            {
                handler->SendSysMessage(LANG_PLAYER_NOT_FOUND);
                handler->SetSentErrorMessage(true);
                return false;
            }
        }
        else
        {
            object = handler->getSelectedUnit();

            if (!object)
            {
                handler->SendSysMessage(LANG_SELECT_CHAR_OR_CREATURE);
                handler->SetSentErrorMessage(true);
                return false;
            }
        }

        CellCoord cellCoord = Trinity::ComputeCellCoord(object->GetPositionX(), object->GetPositionY());
        Cell cell(cellCoord);

        uint32 zoneId, areaId;
        object->GetZoneAndAreaId(zoneId, areaId);

        MapEntry const* mapEntry = sMapStore.LookupEntry(object->GetMapId());
        AreaTableEntry const* zoneEntry = GetAreaEntryByAreaID(zoneId);
        AreaTableEntry const* areaEntry = GetAreaEntryByAreaID(areaId);

        float zoneX = object->GetPositionX();
        float zoneY = object->GetPositionY();

        Map2ZoneCoordinates(zoneX, zoneY, zoneId);

        Map const* map = object->GetMap();
        float groundZ = map->GetHeight(object->GetPhaseMask(), object->GetPositionX(), object->GetPositionY(), MAX_HEIGHT);
        float floorZ = map->GetHeight(object->GetPhaseMask(), object->GetPositionX(), object->GetPositionY(), object->GetPositionZ());

        GridCoord gridCoord = Trinity::ComputeGridCoord(object->GetPositionX(), object->GetPositionY());

        // 63? WHY?
        int gridX = 63 - gridCoord.x_coord;
        int gridY = 63 - gridCoord.y_coord;

        uint32 haveMap = Map::ExistMap(object->GetMapId(), gridX, gridY) ? 1 : 0;
        uint32 haveVMap = Map::ExistVMap(object->GetMapId(), gridX, gridY) ? 1 : 0;

        if (haveVMap)
        {
            if (map->IsOutdoors(object->GetPositionX(), object->GetPositionY(), object->GetPositionZ()))
                handler->PSendSysMessage("You are outdoors");
            else
                handler->PSendSysMessage("You are indoors");
        }
        else
            handler->PSendSysMessage("no VMAP available for area info");

        handler->PSendSysMessage(LANG_MAP_POSITION,
            object->GetMapId(), (mapEntry ? mapEntry->name : "<unknown>"),
            zoneId, (zoneEntry ? zoneEntry->area_name : "<unknown>"),
            areaId, (areaEntry ? areaEntry->area_name : "<unknown>"),
            object->GetPhaseMask(),
            object->GetPositionX(), object->GetPositionY(), object->GetPositionZ(), object->GetOrientation(),
            cell.GridX(), cell.GridY(), cell.CellX(), cell.CellY(), object->GetInstanceId(),
            zoneX, zoneY, groundZ, floorZ, haveMap, haveVMap);

        LiquidData liquidStatus;
        ZLiquidStatus status = map->getLiquidStatus(object->GetPositionX(), object->GetPositionY(), object->GetPositionZ(), MAP_ALL_LIQUIDS, &liquidStatus);

        if (status)
            handler->PSendSysMessage(LANG_LIQUID_STATUS, liquidStatus.level, liquidStatus.depth_level, liquidStatus.type, status);

        return true;
    }

void WorldSession::HandleTextEmoteOpcode(WorldPacket& recvData)
{
    if (!GetPlayer()->IsAlive())
        return;

    if (!GetPlayer()->CanSpeak())
    {
        std::string timeStr = secsToTimeString(m_muteTime - time(NULL));
        SendNotification(GetTrinityString(LANG_WAIT_BEFORE_SPEAKING), timeStr.c_str());
        return;
    }

    uint32 text_emote, emoteNum;
    ObjectGuid guid;

    recvData >> text_emote;
    recvData >> emoteNum;
    recvData >> guid;

    sScriptMgr->OnPlayerTextEmote(GetPlayer(), text_emote, emoteNum, guid);

    EmotesTextEntry const* em = sEmotesTextStore.LookupEntry(text_emote);
    if (!em)
        return;

    uint32 emote_anim = em->textid;

    switch (emote_anim)
    {
        case EMOTE_STATE_SLEEP:
        case EMOTE_STATE_SIT:
        case EMOTE_STATE_KNEEL:
        case EMOTE_ONESHOT_NONE:
            break;
        default:
            // Only allow text-emotes for "dead" entities (feign death included)
            if (GetPlayer()->HasUnitState(UNIT_STATE_DIED))
                break;
            GetPlayer()->HandleEmoteCommand(emote_anim);
            break;
    }

    Unit* unit = ObjectAccessor::GetUnit(*_player, guid);

    CellCoord p = Trinity::ComputeCellCoord(GetPlayer()->GetPositionX(), GetPlayer()->GetPositionY());

    Cell cell(p);
    cell.SetNoCreate();

    Trinity::EmoteChatBuilder emote_builder(*GetPlayer(), text_emote, emoteNum, unit);
    Trinity::LocalizedPacketDo<Trinity::EmoteChatBuilder > emote_do(emote_builder);
    Trinity::PlayerDistWorker<Trinity::LocalizedPacketDo<Trinity::EmoteChatBuilder > > emote_worker(GetPlayer(), sWorld->getFloatConfig(CONFIG_LISTEN_RANGE_TEXTEMOTE), emote_do);
    TypeContainerVisitor<Trinity::PlayerDistWorker<Trinity::LocalizedPacketDo<Trinity::EmoteChatBuilder> >, WorldTypeMapContainer> message(emote_worker);
    cell.Visit(p, message, *GetPlayer()->GetMap(), *GetPlayer(), sWorld->getFloatConfig(CONFIG_LISTEN_RANGE_TEXTEMOTE));

    GetPlayer()->UpdateAchievementCriteria(ACHIEVEMENT_CRITERIA_TYPE_DO_EMOTE, text_emote, 0, unit);

    //Send scripted event call
    if (unit && unit->GetTypeId() == TYPEID_UNIT && ((Creature*)unit)->AI())
        ((Creature*)unit)->AI()->ReceiveEmote(GetPlayer(), text_emote);
}

int immediate_outedges_list(CTXTdeclc callnodeptr call1){
 
  VariantSF subgoal;
  TIFptr tif;
  int j, count = 0,arity;
  Psc psc;
  CPtr oldhreg = NULL;
  struct hashtable *h;	
  struct hashtable_itr *itr;
  callnodeptr cn;
    
  reg[4] = makelist(hreg);
  new_heap_free(hreg);
  new_heap_free(hreg);
  
  if(IsNonNULL(call1)){ /* This can be called from some non incremental predicate */
    h=call1->outedges->hasht;
    
    itr = hashtable1_iterator(h);       
    if (hashtable1_count(h) > 0){
      do {
	cn = hashtable1_iterator_value(itr);
	if(IsNonNULL(cn->goal)){
	  count++;
	  subgoal = (VariantSF) cn->goal;      
	  tif = (TIFptr) subgoal->tif_ptr;
	  psc = TIF_PSC(tif);
	  arity = get_arity(psc);
	  check_glstack_overflow(4,pcreg,2+arity*200); // don't know how much for build_subgoal_args...
	  oldhreg=hreg-2;
	  if(arity>0){
	    sreg = hreg;
	    follow(oldhreg++) = makecs(sreg);
	    hreg += arity + 1;
	    new_heap_functor(sreg, psc);
	    for (j = 1; j <= arity; j++) {
	      new_heap_free(sreg);
	      cell_array1[arity-j] = cell(sreg-1);
	    }
	    build_subgoal_args(subgoal);		
	  }else{
	    follow(oldhreg++) = makestring(get_name(psc));
	  }
	  follow(oldhreg) = makelist(hreg);
	  new_heap_free(hreg);
	  new_heap_free(hreg);
	}
      } while (hashtable1_iterator_advance(itr));
    }
    if (count>0)
      follow(oldhreg) = makenil;
    else
      reg[4] = makenil;
  }else{
    xsb_warn("Called with non-incremental predicate\n");
    reg[4] = makenil;
  }

  //  printterm(stdout,call_list,100);
  return unify(CTXTc reg_term(CTXTc 3),reg_term(CTXTc 4));
}

    void HandleStormSequence(Unit *Cloud) // 1: begin, 2-9: tick, 10: end
    {
        if(StormCount < 10 && StormCount > 1)
        {
            // deal damage
            int32 bp0 = 800;
            for(uint8 i = 2; i < StormCount; ++i)
                bp0 *= 2;

            CellPair p(MaNGOS::ComputeCellPair(m_creature->GetPositionX(), m_creature->GetPositionY()));
            Cell cell(p);
            cell.data.Part.reserved = ALL_DISTRICT;
            cell.SetNoCreate();

            std::list<Unit *> tempUnitMap;

            {
                MaNGOS::AnyAoETargetUnitInObjectRangeCheck u_check(m_creature, m_creature, 999);
                MaNGOS::UnitListSearcher<MaNGOS::AnyAoETargetUnitInObjectRangeCheck> searcher(tempUnitMap, u_check);

                TypeContainerVisitor<MaNGOS::UnitListSearcher<MaNGOS::AnyAoETargetUnitInObjectRangeCheck>, WorldTypeMapContainer > world_unit_searcher(searcher);
                TypeContainerVisitor<MaNGOS::UnitListSearcher<MaNGOS::AnyAoETargetUnitInObjectRangeCheck>, GridTypeMapContainer >  grid_unit_searcher(searcher);

                CellLock<GridReadGuard> cell_lock(cell, p);
                cell_lock->Visit(cell_lock, world_unit_searcher, *(m_creature->GetMap()));
                cell_lock->Visit(cell_lock, grid_unit_searcher, *(m_creature->GetMap()));
            }

            for(std::list<Unit*>::iterator i = tempUnitMap.begin(); i != tempUnitMap.end(); ++i)
            {
                if(!Cloud->IsWithinDistInMap(*i, 15))
                {
                    float x, y, z;
                    (*i)->GetPosition(x, y, z);
                    x = rand()%2 ? x + rand()%5 : x - rand()%5;
                    y = rand()%2 ? y + rand()%5 : y - rand()%5;
                    z = Cloud->GetPositionZ() + 2 - rand()%4; 
                    if(Creature *trigger = m_creature->SummonCreature(MOB_TEMP_TRIGGER, x, y, z, 0, TEMPSUMMON_TIMED_DESPAWN, 2000))
                    {
                        trigger->AddMonsterMoveFlag(MONSTER_MOVE_LEVITATING);
                        trigger->StopMoving();
                        trigger->CastSpell(trigger, 37248, true);
                        trigger->CastCustomSpell(*i, 43137, &bp0, NULL, NULL, true, 0, 0, m_creature->GetGUID());
                    }
                }
            }

            // visual
            float x, y, z;
            for(uint8 i = 0; i < StormCount; ++i)
            {
                Cloud->GetPosition(x, y, z);
                x = rand()%2 ? x + rand()%10 : x - rand()%10;
                y = rand()%2 ? y + rand()%10 : y - rand()%10;
                z = z + 2 - rand()%4; 
                if(Creature *trigger = m_creature->SummonCreature(MOB_TEMP_TRIGGER, x, y, z, 0, TEMPSUMMON_TIMED_DESPAWN, 2000))
                {
                    trigger->AddMonsterMoveFlag(MONSTER_MOVE_LEVITATING);
                    trigger->StopMoving();
                    trigger->CastSpell(trigger, 37248, true);
                }
                Cloud->GetPosition(x, y, z);
                x = rand()%2 ? x + 10 + rand()%10 : x - 10 - rand()%10;
                y = rand()%2 ? y + 10 + rand()%10 : y - 10 - rand()%10;
                if(Creature *trigger = m_creature->SummonCreature(MOB_TEMP_TRIGGER, x, y, m_creature->GetPositionZ(), 0, TEMPSUMMON_TIMED_DESPAWN, 2000))
                {
                    trigger->SetMaxHealth(9999999);
                    trigger->SetHealth(9999999);
                    trigger->CastSpell(trigger, 43661, true);
                }
            }
        }

        StormCount++;
        if(StormCount > 10)
        {
            StormCount = 0; // finish
            SummonEagles_Timer = 5000;
            m_creature->InterruptNonMeleeSpells(false);
            Cloud->RemoveAurasDueToSpell(45213);
            CloudGUID = 0;
            if(Unit* Cyclone = Unit::GetUnit(*m_creature, CycloneGUID))
                Cyclone->RemoveAurasDueToSpell(25160);
            SetWeather(WEATHER_STATE_FINE, 0.0f);
            isRaining = false;
        }

        StormSequenceTimer = 1000;
    }

int main(int argc, char* argv[])
{
  int i, j, loop, num_alive, maxloop;
  int ldboard, ldnbngb, ldlboard;
  double t1, t2;
  double temps;
 
  int *board;
  int *nbngb;

  int local_alive;
  int *global_board;

  struct grid grid;
  MPI_Comm comm;
  int nb_proc_row;
  int nb_proc_tot;
  int rank;
  int nb_in_block;

  MPI_Init(&argc,&argv);
  MPI_Comm_size(MPI_COMM_WORLD, &nb_proc_tot);
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);

  // initialization of the grid communicator
  if (EXIT_FAILURE == compute_communicator(nb_proc_tot,&nb_proc_row,&comm,&rank)){
    MPI_Finalize();
    return EXIT_SUCCESS;
  }


  if (argc < 2) {
    maxloop = 10;
  } else if (argc > 2){
    maxloop = atoi(argv[1]);
    BS = atoi(argv[2]);
  } else
    maxloop = atoi(argv[1]);
  num_alive = 0;
  local_alive = 0;

  /* Leading dimension of the board array */
  ldboard = BS;
  if (ldboard % nb_proc_row != 0){
    if (rank == 0)
      printf("Wrong BS (or wrong number of procs) ... exiting now.\n");
    MPI_Finalize();
    return EXIT_FAILURE;
  }

  /* Leading dimension of the neigbour counters array */
  nb_in_block = ldboard / nb_proc_row;
  ldnbngb = nb_in_block;
  ldlboard = nb_in_block + 2;

  board = malloc( ldlboard * ldlboard * sizeof(int) );
  nbngb = malloc( ldnbngb * ldnbngb * sizeof(int) );

  if (rank == 0){
    global_board = malloc( ldboard * ldboard * sizeof(int) );
    num_alive = generate_initial_board( &global_cell( 1, 1), ldboard );
    printf("Starting number of living cells = %d\n", num_alive);
    t1 = mytimer();
  }

  matrix_placement_proc(nb_proc_row, nb_in_block, &comm, &(global_cell( 1, 1)), &(cell( 1, 1)), SCATTER, ldlboard);

  mpi_grid_init(&comm, &grid, rank);
  //printf("rank #%d: %d %d\n", rank, grid.rank_I, grid.rank_J);


  //output_lboard( nb_in_block, board, ldlboard, 0, rank );

  for (loop = 1; loop <= maxloop; loop++) {

    MPI_Datatype blocktype; // we need a specific type for row exchange
    MPI_Type_vector(nb_in_block, 1, ldlboard, MPI_INT, &blocktype);
    MPI_Type_commit(&blocktype);
    // for upper/lower ghost row
    MPI_Sendrecv(&(cell( 1, 1)), 1, blocktype, grid.proc_above, 99, 
		 &(cell( nb_in_block+1, 1)), 1, blocktype, grid.proc_under, 99,
		 comm, MPI_STATUS_IGNORE);
    MPI_Sendrecv(&(cell( nb_in_block, 1)), 1, blocktype, grid.proc_under, 99,
		 &(cell( 0, 1)), 1, blocktype, grid.proc_above, 99, 
		 comm, MPI_STATUS_IGNORE);

    // for left/right ghost col
    MPI_Sendrecv(&(cell( 0, 1)), ldlboard, MPI_INT, grid.proc_left, 98, 
		 &(cell( 0, nb_in_block+1)), ldlboard, MPI_INT, grid.proc_right, 98,
		 comm, MPI_STATUS_IGNORE);
    MPI_Sendrecv(&(cell( 0, nb_in_block)), ldlboard, MPI_INT, grid.proc_right, 98,
		 &(cell( 0, 0)), ldlboard, MPI_INT, grid.proc_left, 98, 
		 comm, MPI_STATUS_IGNORE);

    //debug
    /* if (loop == 1) */
    /*   output_lboard( nb_in_block, board, ldlboard, 0, rank ); */

    //calcul du nombre de voisins
    for (j = 1; j <= nb_in_block; j++) {
//.........这里部分代码省略.........

TPointD XshHandleManager::getHandlePos(const TStageObjectId &id,
                                       const std::string &handle,
                                       int row) const {
  static const double unit = 8.0;

  assert(m_xsh->getScene());

  if (handle == "")
    return TPointD();

  else if (handle[0] == 'H' && handle.length() > 1) {
    // Hook port case

    if (!id.isColumn()) return TPointD();

    int hookIndex = 0;
    {
      int h, hLength = handle.length();
      for (h = 1; h != hLength; ++h)
        hookIndex = hookIndex * 10 + (handle[h] - '0');
    }

    TStageObject *obj = m_xsh->getStageObject(id);
    if (const PlasticSkeletonDeformationP &def =
            obj->getPlasticSkeletonDeformation()) {
      int skelId = def->skeletonId(row);

      PlasticSkeleton skel;
      def->storeDeformedSkeleton(skelId, row, skel);

      int v = def->vertexIndex(hookIndex, skelId);
      return (v >= 0) ? TScale(1.0 / Stage::inch) * skel.vertex(v).P()
                      : TPointD();
    }

    --hookIndex;

    int col = id.getIndex();
    TXshCell cell(m_xsh->getCell(row, col));

    TXshLevel *xl = cell.m_level.getPointer();
    if (!xl) return TPointD();

    TXshSimpleLevel *sl = xl->getSimpleLevel();
    if (!sl) return TPointD();

    Hook *hook = sl->getHookSet()->getHook(hookIndex);
    if (!hook) return TPointD();

    TFrameId fid(cell.m_frameId);

    TPointD pos = hook->getAPos(fid) * (1.0 / Stage::inch);
    TPointD delta;

    for (int r = row - 1; r >= 0; --r) {
      cell = m_xsh->getCell(r, col);
      if (cell.m_level.getPointer() != xl) break;

      const TFrameId &precFid = cell.m_frameId;
      delta += computePassHook(fid, precFid, sl, hook);
      fid = precFid;
    }

    pos += delta * (1.0 / Stage::inch);
    pos = getDpiAffine(sl, cell.m_frameId, true) * pos;

    return pos;
  } else if ('A' <= handle[0] && handle[0] <= 'Z')
    return TPointD(unit * (handle[0] - 'B'), 0);
  else if ('a' <= handle[0] && handle[0] <= 'z')
    return TPointD(0.5 * unit * (handle[0] - 'b'), 0);
  else
    return TPointD();
}

const cell worksheet::get_cell(const cell_reference &reference) const
{
    return cell(&d_->cell_map_.at(reference.get_row()).at(reference.get_column_index()));
}

void
tbl_layout(struct tbl_node *tbl, int ln, const char *p, int pos)
{
	struct tbl_row	*rp;

	rp = NULL;
	for (;;) {
		/* Skip whitespace before and after each cell. */

		while (p[pos] == ' ' || p[pos] == '\t')
			pos++;

		switch (p[pos]) {
		case ',':  /* Next row on this input line. */
			pos++;
			rp = NULL;
			continue;
		case '\0':  /* Next row on next input line. */
			return;
		case '.':  /* End of layout. */
			pos++;
			tbl->part = TBL_PART_DATA;

			/*
			 * When the layout is completely empty,
			 * default to one left-justified column.
			 */

			if (tbl->first_row == NULL) {
				tbl->first_row = tbl->last_row =
				    mandoc_calloc(1, sizeof(*rp));
			}
			if (tbl->first_row->first == NULL) {
				mandoc_msg(MANDOCERR_TBLLAYOUT_NONE,
				    tbl->parse, ln, pos, NULL);
				cell_alloc(tbl, tbl->first_row,
				    TBL_CELL_LEFT);
				return;
			}

			/*
			 * Search for the widest line
			 * along the left and right margins.
			 */

			for (rp = tbl->first_row; rp; rp = rp->next) {
				if (tbl->opts.lvert < rp->vert)
					tbl->opts.lvert = rp->vert;
				if (rp->last != NULL &&
				    rp->last->col + 1 == tbl->opts.cols &&
				    tbl->opts.rvert < rp->last->vert)
					tbl->opts.rvert = rp->last->vert;

				/* If the last line is empty, drop it. */

				if (rp->next != NULL &&
				    rp->next->first == NULL) {
					free(rp->next);
					rp->next = NULL;
				}
			}
			return;
		default:  /* Cell. */
			break;
		}

		/*
		 * If the last line had at least one cell,
		 * start a new one; otherwise, continue it.
		 */

		if (rp == NULL) {
			if (tbl->last_row == NULL ||
			    tbl->last_row->first != NULL) {
				rp = mandoc_calloc(1, sizeof(*rp));
				if (tbl->last_row)
					tbl->last_row->next = rp;
				else
					tbl->first_row = rp;
				tbl->last_row = rp;
			} else
				rp = tbl->last_row;
		}
		cell(tbl, rp, ln, p, &pos);
	}
}

template<> void ObjectGridLoader::SetObjectCell(Creature* obj, CellPair const& cell_pair)
{
    Cell cell(cell_pair);

    obj->SetCurrentCell(cell);
}

int main() {
	double dx = 1.0 / 100;
    double Pinf_water = 1.1645*1000000000;
    double gama_water = 1.932;

	std::vector<CSE::Cell> cell(2);
     
	std::map<std::string, double> leftCell_phase1, leftCell_phase2, rightCell_phase1, rightCell_phase2;
	leftCell_phase1["gama"] = 1.4; 
	leftCell_phase1["density"] = 112.99;
	leftCell_phase1["pressure"] = 1e+07;
	leftCell_phase1["velocity"] = 0;
	leftCell_phase1["phi"] = 0.0;
    leftCell_phase2["Pinf"] = 0;

	leftCell_phase2["gama"] = gama_water;
	leftCell_phase2["density"] = 1025;
	leftCell_phase2["pressure"] = 1e+07;
	leftCell_phase2["velocity"] = 0;
	leftCell_phase2["phi"] = 1-leftCell_phase1["phi"];
    leftCell_phase2["Pinf"] = Pinf_water;

	cell[0].addPhase(leftCell_phase1);
	cell[0].addPhase(leftCell_phase2);

	rightCell_phase1["gama"] = 1.4;
	rightCell_phase1["density"] = 56.49;
	rightCell_phase1["pressure"] = 5e+06;
	rightCell_phase1["velocity"] = 0;
	rightCell_phase1["phi"] = 1.0;
    rightCell_phase2["Pinf"] = 0;

	rightCell_phase2["gama"] = gama_water;
	rightCell_phase2["density"] = 1025;
	rightCell_phase2["pressure"] = 5e+06;
	rightCell_phase2["velocity"] = 0;
	rightCell_phase2["phi"] = 1-rightCell_phase1["phi"];
    rightCell_phase2["Pinf"] = Pinf_water;

	cell[1].addPhase(rightCell_phase1);
	cell[1].addPhase(rightCell_phase2);
	//for (unsigned i = 0; i < 30; ++i)
	//	cell[i].addPhase(leftCell_phase1);
	//for (unsigned i = 30; i < 102; ++i)
	//	cell[i].addPhase(rightCell_phase1);
//	for (unsigned i = 0; i < 50; ++i)
//			cell[i].addPhase(1.4, 1, 0.4, -2.0, 0);
//    for (unsigned i = 50; i < 102; ++i)
//			cell[i].addPhase(1.4, 1, 0.4, 2.0, 0);
	for (unsigned i = 0; i < cell.size(); ++i)
		cell[i].computeConsVar();
//
//	CSE::FluxComputation *fluxMethod = new CSE::Godunov(cell, 0.2);
//	fluxMethod->Solve(0.2);
//	fluxMethod->copyResult(cell);
//	std::string fluxName;
//
//	fluxName = variableType<CSE::FluxComputation>(*fluxMethod);
//	CSE::printPhaseValues(fluxName, cell);

    CSE::RiemannSolver solver(cell[0], cell[1], 2, false);
	solver.multiphaseSolve();
    std::cout << "Over" << std::endl;
	return 1;
}

uint32_t MxField::monteCarlo(uint32_t stepMax) const
{
    uint32_t positive = 0;
    
    double rmin = NitroDiameter / 2.0;
//    std::vector<Cell> clusters = cells();
//    std::cout << clusters.size() << std::endl;

//    std::vector<std::vector<std::vector<std::vector<Cell> > > > grid;
//    grid.resize(m_sides.x);
//    for (uint32_t x = 0; x < m_sides.x; ++x) {
//        grid[x].resize(m_sides.y);
//        for (uint32_t y = 0; y < m_sides.y; ++y) {
//            grid[x][y].resize(m_sides.z);
//        }
//    }
//    for (const Cell& cell : clusters) {
//        grid[uint32_t(cell.coord().x * m_sides.x / m_sizes.x)]
//                [uint32_t(cell.coord().y * m_sides.y / m_sizes.y)]
//                [uint32_t(cell.coord().z * m_sides.z / m_sizes.z)].push_back(cell);
//    }

    for (uint32_t i = 0; i < stepMax; ) {
        int32_t rx = rand() % int32_t(m_sides.x);
        int32_t ry = rand() % int32_t(m_sides.y);
        int32_t rz = rand() % int32_t(m_sides.z);
        if (m_field[size_t(rx)][size_t(ry)][size_t(rz)] != 0) {
            Cell curr(new FSphere(radius()), dCoord(rx * side(), ry * side(), rz * side()));
//            uint32_t rcluster = rand() % (uint32_t(clusters.size()));
//            const Cell& curr = clusters[rcluster];

            double xc = curr.coord().x;
            double yc = curr.coord().y;
            double zc = curr.coord().z;
            double rc = curr.figure()->radius();

            //spheric!
            double teta = 2.0 * M_PI * (rand() / double(RAND_MAX));
            double phi  = 2.0 * M_PI * (rand() / double(RAND_MAX));

            double ixc = xc + (rc + rmin) * sin(teta) * cos(phi);
            double iyc = yc + (rc + rmin) * sin(teta) * sin(phi);
            double izc = zc + (rc + rmin) * cos(teta);

            Cell cell(new FSphere(rmin), dCoord(ixc, iyc, izc));

            bool overlap = false;
            for (int32_t ix = rx - 2; ix < rx + 2; ++ix) {
                if (overlap) {
                    break;
                }
                for (int32_t iy = ry - 2; iy < ry + 2; ++iy) {
                    if (overlap) {
                        break;
                    }
                    for (int32_t iz = rz - 2; iz < rz + 2; ++iz) {
                        if (ix != rx && iy != ry && iz != rz) {
                            Cell temp(new FSphere(radius()),
                                      dCoord((uint32_t(ix + int32_t(m_sides.x)) % m_sides.x) * side(),
                                             (uint32_t(iy + int32_t(m_sides.y)) % m_sides.y) * side(),
                                             (uint32_t(iz + int32_t(m_sides.z)) % m_sides.z) * side()));
                            if (isOverlapped(&temp, &cell)) {
                                overlap = true;
                                break;
                            }
                        }
                    }
                }
            }
//            for (size_t ic = 0; ic < clusters.size(); ++ic) {
//                if (ic != rcluster) {
//                    if (is_overlapped(clusters[ic], cell)) {
//                        overlap = true;
//                        break;
//                    }
//                }
//            }
            if (!overlap) {
                ++positive;
            }
            ++i;
        }
    }
    return positive;
}

cell proc_mul(const cells & c)
{
    long n(1);
    for (cellit i = c.begin(); i != c.end(); ++i) n *= atol(i->val.c_str());
    return cell(Number, str(n));
}

int create_lazy_call_list(CTXTdeclc  callnodeptr call1){
  VariantSF subgoal;
  TIFptr tif;
  int j,count=0,arity; 
  Psc psc;
  CPtr oldhreg=NULL;

  //  print_call_list(lazy_affected);

  reg[6] = reg[5] = makelist(hreg);  // reg 5 first not-used, use regs in case of stack expanson
  new_heap_free(hreg);   // make heap consistent
  new_heap_free(hreg);
  while((call1 = delete_calllist_elt(&lazy_affected)) != EMPTY){
    subgoal = (VariantSF) call1->goal;      
    //    fprintf(stddbg,"  considering ");print_subgoal(stdout,subgoal);printf("\n");
    if(IsNULL(subgoal)){ /* fact predicates */
      call1->deleted = 0; 
      continue;
    }
    if (subg_visitors(subgoal)) {
      sprint_subgoal(CTXTc forest_log_buffer_1,0,subgoal);
#ifdef WARN_ON_UNSAFE_UPDATE
      xsb_warn("%d Choice point(s) exist to the table for %s -- cannot incrementally update (create_lazy_call_list)\n",
	       subg_visitors(subgoal),forest_log_buffer_1->fl_buffer);
#else
      xsb_abort("%d Choice point(s) exist to the table for %s -- cannot incrementally update (create_lazy_call_list)\n",
	       subg_visitors(subgoal),forest_log_buffer_1->fl_buffer);
#endif
      continue;
    }
    //    fprintf(stddbg,"adding dependency for ");print_subgoal(stdout,subgoal);printf("\n");

    count++;
    tif = (TIFptr) subgoal->tif_ptr;
    //    if (!(psc = TIF_PSC(tif)))
    //	xsb_table_error(CTXTc "Cannot access dynamic incremental table\n");	
    psc = TIF_PSC(tif);
    arity = get_arity(psc);
    check_glstack_overflow(6,pcreg,2+arity*200); // don't know how much for build_subgoal_args...
    oldhreg = clref_val(reg[6]);  // maybe updated by re-alloc
    if(arity>0){
      sreg = hreg;
      follow(oldhreg++) = makecs(sreg);
      hreg += arity + 1;  // had 10, why 10?  why not 3? 2 for list, 1 for functor (dsw)
      new_heap_functor(sreg, psc);
      for (j = 1; j <= arity; j++) {
	new_heap_free(sreg);
	cell_array1[arity-j] = cell(sreg-1);
      }
      build_subgoal_args(subgoal);		
    } else {
      follow(oldhreg++) = makestring(get_name(psc));
    }
    reg[6] = follow(oldhreg) = makelist(hreg);
    new_heap_free(hreg);
    new_heap_free(hreg);
  }
  if(count > 0) {
    follow(oldhreg) = makenil;
    hreg -= 2;  /* take back the extra words allocated... */
  } else
    reg[5] = makenil;
    
  return unify(CTXTc reg_term(CTXTc 4),reg_term(CTXTc 5));

  /*int i;
    for(i=0;i<callqptr;i++){
      if(IsNonNULL(callq[i]) && (callq[i]->deleted==1)){
    sfPrintGoal(stdout,(VariantSF)callq[i]->goal,NO);
    printf(" %d %d\n",callq[i]->falsecount,callq[i]->deleted);
    }
    }
  printf("-----------------------------\n");   */
}

cell proc_length(const cells & c) { return cell(Number, str(c[0].list.size())); }

// define the bare minimum set of primintives necessary to pass the unit tests
void add_globals(environment & env)
{
    env["nil"] = nil;   env["#f"] = false_sym;  env["#t"] = true_sym;
    env["append"] = cell(&proc_append);   env["car"]  = cell(&proc_car);
    env["cdr"]    = cell(&proc_cdr);      env["cons"] = cell(&proc_cons);
    env["length"] = cell(&proc_length);   env["list"] = cell(&proc_list);
    env["null?"]  = cell(&proc_nullp);    env["+"]    = cell(&proc_add);
    env["-"]      = cell(&proc_sub);      env["*"]    = cell(&proc_mul);
    env["/"]      = cell(&proc_div);      env[">"]    = cell(&proc_greater);
    env["<"]      = cell(&proc_less);     env["<="]   = cell(&proc_less_equal);
}

void mitk::ReduceContourSetFilter::GenerateData()
{
  unsigned int numberOfInputs = this->GetNumberOfInputs();
  unsigned int numberOfOutputs (0);

  vtkSmartPointer<vtkPolyData> newPolyData;
  vtkSmartPointer<vtkCellArray> newPolygons;
  vtkSmartPointer<vtkPoints> newPoints;

  //For the purpose of evaluation
//  unsigned int numberOfPointsBefore (0);
  m_NumberOfPointsAfterReduction=0;

  for(unsigned int i = 0; i < numberOfInputs; i++)
  {
    mitk::Surface* currentSurface = const_cast<mitk::Surface*>( this->GetInput(i) );
    vtkSmartPointer<vtkPolyData> polyData = currentSurface->GetVtkPolyData();

    newPolyData = vtkPolyData::New();
    newPolygons = vtkCellArray::New();
    newPoints = vtkPoints::New();

    vtkSmartPointer<vtkCellArray> existingPolys = polyData->GetPolys();

    vtkSmartPointer<vtkPoints> existingPoints = polyData->GetPoints();

    existingPolys->InitTraversal();

    vtkIdType* cell (NULL);
    vtkIdType cellSize (0);

    for( existingPolys->InitTraversal(); existingPolys->GetNextCell(cellSize, ce