C++ Cells类代码示例

    /// Gas viscosity.
    /// \param[in]  pg     Array of n gas pressure values.
    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
    /// \return            Array of n viscosity values.
    ADB BlackoilPropsAd::muGas(const ADB& pg,
                               const Cells& cells) const
    {
#if 1
        return ADB::constant(muGas(pg.value(), cells), pg.blockPattern());
#else
        if (!pu_.phase_used[Gas]) {
            THROW("Cannot call muGas(): gas phase not present.");
        }
        const int n = cells.size();
        ASSERT(pg.value().size() == n);
        const int np = props_.numPhases();
        Block z = Block::Zero(n, np);
        Block mu(n, np);
        Block dmu(n, np);
        props_.viscosity(n, pg.value().data(), z.data(), cells.data(), mu.data(), dmu.data());
        ADB::M dmu_diag = spdiag(dmu.col(pu_.phase_pos[Gas]));
        const int num_blocks = pg.numBlocks();
        std::vector<ADB::M> jacs(num_blocks);
        for (int block = 0; block < num_blocks; ++block) {
            jacs[block] = dmu_diag * pg.derivative()[block];
        }
        return ADB::function(mu.col(pu_.phase_pos[Gas]), jacs);
#endif
    }

void GuiPart::updateConnectorPoints( bool add )
{
	ICNDocument *icnd = dynamic_cast<ICNDocument*>(p_parent->itemDocument());
	if ( !icnd)
		return;
	
	Cells * cells = icnd->cells();
	if (!cells)
		return;
	
	if ( !isVisible() )
		add = false;
	
	if ( add == b_pointsAdded )
		return;
	
	b_pointsAdded = add;
	
	int mult = add ? 1 : -1;
	int sx = roundDown( x(), 8 );
	int sy = roundDown( y(), 8 );
	int ex = roundDown( x()+width(), 8 );
	int ey = roundDown( y()+height(), 8 );
	
	for ( int x=sx; x<=ex; ++x )
	{
		for ( int y=sy; y<=ey; ++y )
		{
			if ( cells->haveCell( x, y ) )
				cells->cell( x, y ).CIpenalty += mult*ICNDocument::hs_item/2;
		}
	}
}

// The strategy here is, though this is n^2, it's a very small n, just
// finely-spaced adjacent cells.
bool SpatialHash::isWithinDistanceOfAnything(PointType const& physPt,
                                             double distance) const
{
   // (save taking a lot of square roots)
   double d2 = distance * distance;

   Index idx = index_of(physPt);
   Cells nbrs;
   get_neighbors(idx, nbrs);
   // 27 <= we include the center cell itself, too.
   assert(nbrs.size() <= 27);
   for (Cells::const_iterator itCells = nbrs.begin(), endCells = nbrs.end();
        itCells != endCells; ++itCells) {
      Pts const* pts = *itCells;
      if (pts) {
         for (Pts::const_iterator itPts = pts->begin(), endPts = pts->end();
              itPts != endPts; ++itPts) {
            if (itPts->SquaredEuclideanDistanceTo<double>(physPt) < d2) {
               return true;
            }
         }
      }
   }
   return false;
}

    /// Water viscosity.
    /// \param[in]  pw     Array of n water pressure values.
    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
    /// \return            Array of n viscosity values.
    ADB BlackoilPropsAd::muWat(const ADB& pw,
                               const Cells& cells) const
    {
#if 1
        return ADB::constant(muWat(pw.value(), cells), pw.blockPattern());
#else
        if (!pu_.phase_used[Water]) {
            OPM_THROW(std::runtime_error, "Cannot call muWat(): water phase not present.");
        }
        const int n = cells.size();
        assert(pw.value().size() == n);
        const int np = props_.numPhases();
        Block z = Block::Zero(n, np);
        Block mu(n, np);
        Block dmu(n, np);
        props_.viscosity(n, pw.value().data(), z.data(), cells.data(), mu.data(), dmu.data());
        ADB::M dmu_diag = spdiag(dmu.col(pu_.phase_pos[Water]));
        const int num_blocks = pw.numBlocks();
        std::vector<ADB::M> jacs(num_blocks);
        for (int block = 0; block < num_blocks; ++block) {
            jacs[block] = dmu_diag * pw.derivative()[block];
        }
        return ADB::function(mu.col(pu_.phase_pos[Water]), jacs);
#endif
    }

void VTKWriter::initializeOutput(int numParticles) {

	vtkFile = new VTKFile_t("UnstructuredGrid");

	// per point, we add type, position, velocity and force
	PointData pointData;
	DataArray_t mass(type::Float32, "mass", 1);
	DataArray_t velocity(type::Float32, "velocity", 3);
	DataArray_t forces(type::Float32, "force", 3);
	DataArray_t type(type::Int32, "type", 1);
	pointData.DataArray().push_back(mass);
	pointData.DataArray().push_back(velocity);
    pointData.DataArray().push_back(forces);
    pointData.DataArray().push_back(type);

	CellData cellData; // we don't have cell data => leave it empty

	// 3 coordinates
	Points points;
	DataArray_t pointCoordinates(type::Float32, "points", 3);
	points.DataArray().push_back(pointCoordinates);

	Cells cells; // we don't have cells, => leave it empty
	// for some reasons, we have to add a dummy entry for paraview
	DataArray_t cells_data(type::Float32, "types", 0);
	cells.DataArray().push_back(cells_data);

	PieceUnstructuredGrid_t piece(pointData, cellData, points, cells, numParticles, 0);
	UnstructuredGrid_t unstructuredGrid(piece);
	vtkFile->UnstructuredGrid(unstructuredGrid);
}

 /// Oil formation volume factor.
 /// \param[in]  po     Array of n oil pressure values.
 /// \param[in]  rs     Array of n gas solution factor values.
 /// \param[in]  cond   Array of n taxonomies classifying fluid condition.
 /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
 /// \return            Array of n formation volume factor values.
 ADB BlackoilPropsAd::bOil(const ADB& po,
                           const ADB& rs,
                           const std::vector<PhasePresence>& /*cond*/,
                           const Cells& cells) const
 {
     if (!pu_.phase_used[Oil]) {
         OPM_THROW(std::runtime_error, "Cannot call muOil(): oil phase not present.");
     }
     const int n = cells.size();
     assert(po.value().size() == n);
     const int np = props_.numPhases();
     Block z = Block::Zero(n, np);
     if (pu_.phase_used[Gas]) {
         // Faking a z with the right ratio:
         //   rs = zg/zo
         z.col(pu_.phase_pos[Oil]) = V::Ones(n, 1);
         z.col(pu_.phase_pos[Gas]) = rs.value();
     }
     Block matrix(n, np*np);
     Block dmatrix(n, np*np);
     props_.matrix(n, po.value().data(), z.data(), cells.data(), matrix.data(), dmatrix.data());
     const int phase_ind = pu_.phase_pos[Oil];
     const int column = phase_ind*np + phase_ind; // Index of our sought diagonal column.
     ADB::M db_diag = spdiag(dmatrix.col(column));
     const int num_blocks = po.numBlocks();
     std::vector<ADB::M> jacs(num_blocks);
     for (int block = 0; block < num_blocks; ++block) {
         // For now, we deliberately ignore the derivative with respect to rs,
         // since the BlackoilPropertiesInterface class does not evaluate it.
         // We would add to the next line: + db_drs_diag * rs.derivative()[block]
         jacs[block] = db_diag * po.derivative()[block];
     }
     return ADB::function(matrix.col(column), jacs);
 }

    /// Gas viscosity.
    /// \param[in]  pg     Array of n gas pressure values.
    /// \param[in]  rv     Array of n vapor oil/gas ratio
    /// \param[in]  cond   Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
    /// \return            Array of n viscosity values.
    ADB BlackoilPropsAd::muGas(const ADB& pg,
                               const ADB& rv,
                               const std::vector<PhasePresence>& cond,
                               const Cells& cells) const
    {
#if 1
        return ADB::constant(muGas(pg.value(), rv.value(),cond,cells), pg.blockPattern());
#else
        if (!pu_.phase_used[Gas]) {
            OPM_THROW(std::runtime_error, "Cannot call muGas(): gas phase not present.");
        }
        const int n = cells.size();
        assert(pg.value().size() == n);
        const int np = props_.numPhases();
        Block z = Block::Zero(n, np);
        if (pu_.phase_used[Oil]) {
            // Faking a z with the right ratio:
            //   rv = zo/zg
            z.col(pu_.phase_pos[Oil]) = rv;
            z.col(pu_.phase_pos[Gas]) = V::Ones(n, 1);
        }
        Block mu(n, np);
        Block dmu(n, np);
        props_.viscosity(n, pg.value().data(), z.data(), cells.data(), mu.data(), dmu.data());
        ADB::M dmu_diag = spdiag(dmu.col(pu_.phase_pos[Gas]));
        const int num_blocks = pg.numBlocks();
        std::vector<ADB::M> jacs(num_blocks);
        for (int block = 0; block < num_blocks; ++block) {
            jacs[block] = dmu_diag * pg.derivative()[block];
        }
        return ADB::function(mu.col(pu_.phase_pos[Gas]), jacs);
#endif
    }

 /// Gas formation volume factor.
 /// \param[in]  pg     Array of n gas pressure values.
 /// \param[in]  rv     Array of n vapor oil/gas ratio
 /// \param[in]  cond   Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
 /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
 /// \return            Array of n formation volume factor values.
 ADB BlackoilPropsAd::bGas(const ADB& pg,
                           const ADB& rv,
                           const std::vector<PhasePresence>& /*cond*/,
                           const Cells& cells) const
 {
     if (!pu_.phase_used[Gas]) {
         OPM_THROW(std::runtime_error, "Cannot call muGas(): gas phase not present.");
     }
     const int n = cells.size();
     assert(pg.value().size() == n);
     const int np = props_.numPhases();
     Block z = Block::Zero(n, np);
     if (pu_.phase_used[Oil]) {
         // Faking a z with the right ratio:
         //   rv = zo/zg
         z.col(pu_.phase_pos[Oil]) = rv.value();
         z.col(pu_.phase_pos[Gas]) = V::Ones(n, 1);
     }
     Block matrix(n, np*np);
     Block dmatrix(n, np*np);
     props_.matrix(n, pg.value().data(), z.data(), cells.data(), matrix.data(), dmatrix.data());
     const int phase_ind = pu_.phase_pos[Gas];
     const int column = phase_ind*np + phase_ind; // Index of our sought diagonal column.
     ADB::M db_diag = spdiag(dmatrix.col(column));
     const int num_blocks = pg.numBlocks();
     std::vector<ADB::M> jacs(num_blocks);
     for (int block = 0; block < num_blocks; ++block) {
         jacs[block] = db_diag * pg.derivative()[block];
     }
     return ADB::function(matrix.col(column), jacs);
 }

 /// Relative permeabilities for all phases.
 /// \param[in]  sw     Array of n water saturation values.
 /// \param[in]  so     Array of n oil saturation values.
 /// \param[in]  sg     Array of n gas saturation values.
 /// \param[in]  cells  Array of n cell indices to be associated with the saturation values.
 /// \return            An std::vector with 3 elements, each an array of n relperm values,
 ///                    containing krw, kro, krg. Use PhaseIndex for indexing into the result.
 std::vector<V> BlackoilPropsAd::relperm(const V& sw,
                                         const V& so,
                                         const V& sg,
                                         const Cells& cells) const
 {
     const int n = cells.size();
     const int np = props_.numPhases();
     Block s_all(n, np);
     if (pu_.phase_used[Water]) {
         assert(sw.size() == n);
         s_all.col(pu_.phase_pos[Water]) = sw;
     }
     if (pu_.phase_used[Oil]) {
         assert(so.size() == n);
         s_all.col(pu_.phase_pos[Oil]) = so;
     }
     if (pu_.phase_used[Gas]) {
         assert(sg.size() == n);
         s_all.col(pu_.phase_pos[Gas]) = sg;
     }
     Block kr(n, np);
     props_.relperm(n, s_all.data(), cells.data(), kr.data(), 0);
     std::vector<V> relperms;
     relperms.reserve(3);
     for (int phase = 0; phase < 3; ++phase) {
         if (pu_.phase_used[phase]) {
             relperms.emplace_back(kr.col(pu_.phase_pos[phase]));
         } else {
             relperms.emplace_back();
         }
     }
     return relperms;
 }

    /// Oil viscosity.
    /// \param[in]  po     Array of n oil pressure values.
    /// \param[in]  rs     Array of n gas solution factor values.
    /// \param[in]  cells  Array of n cell indices to be associated with the pressure values.
    /// \return            Array of n viscosity values.
    ADB BlackoilPropsAd::muOil(const ADB& po,
                               const ADB& rs,
                               const Cells& cells) const
    {
#if 1
        return ADB::constant(muOil(po.value(), rs.value(), cells), po.blockPattern());
#else
        if (!pu_.phase_used[Oil]) {
            THROW("Cannot call muOil(): oil phase not present.");
        }
        const int n = cells.size();
        ASSERT(po.value().size() == n);
        const int np = props_.numPhases();
        Block z = Block::Zero(n, np);
        if (pu_.phase_used[Gas]) {
            // Faking a z with the right ratio:
            //   rs = zg/zo
            z.col(pu_.phase_pos[Oil]) = V::Ones(n, 1);
            z.col(pu_.phase_pos[Gas]) = rs.value();
        }
        Block mu(n, np);
        Block dmu(n, np);
        props_.viscosity(n, po.value().data(), z.data(), cells.data(), mu.data(), dmu.data());
        ADB::M dmu_diag = spdiag(dmu.col(pu_.phase_pos[Oil]));
        const int num_blocks = po.numBlocks();
        std::vector<ADB::M> jacs(num_blocks);
        for (int block = 0; block < num_blocks; ++block) {
            // For now, we deliberately ignore the derivative with respect to rs,
            // since the BlackoilPropertiesInterface class does not evaluate it.
            // We would add to the next line: + dmu_drs_diag * rs.derivative()[block]
            jacs[block] = dmu_diag * po.derivative()[block];
        }
        return ADB::function(mu.col(pu_.phase_pos[Oil]), jacs);
#endif
    }

/** TODO: there's redundancy in adding the data members (the same has to be done
 * for the parallel VTK File. -> use array of structs?
 */
void VTKGridWriterImplementation::initializeVTKFile() {
	PointData pointData;
	// we don't need point data at all!?
	//DataArray_t position(type::Int32, "id", 0);
	//pointData.DataArray().push_back(position);

	CellData cellData;
	DataArray_t cells_count(type::Int32, "numberOfMolecules", 1);
	cellData.DataArray().push_back(cells_count);
	DataArray_t node_rank(type::Int32, "node-rank", 1);
	cellData.DataArray().push_back(node_rank);
	DataArray_t index(type::UInt32, "index", 1);
	cellData.DataArray().push_back(index);
	DataArray_t velocity(type::Float64, "velocity", 3);
	cellData.DataArray().push_back(velocity);

	// 3 coordinates
	Points points;
	DataArray_t pointCoordinates(type::Float32, "points", 3);
	points.DataArray().push_back(pointCoordinates);

	Cells cells;
	DataArray_t cells_connectivity(type::Int32, "connectivity", 1);
	cells.DataArray().push_back(cells_connectivity);
	DataArray_t cells_offsets(type::Int32, "offsets", 1);
	cells.DataArray().push_back(cells_offsets);
	DataArray_t cells_type(type::Int32, "types", 1);
	cells.DataArray().push_back(cells_type);

	PieceUnstructuredGrid_t piece(pointData, cellData, points, cells, 0, 0);
	UnstructuredGrid_t unstructuredGrid(piece);
	_vtkFile = new VTKFile_t("UnstructuredGrid");
	_vtkFile->UnstructuredGrid(unstructuredGrid);
}

int main(int argc, char *argv[])
{
    Cells cells;

    cells.run();

	return 0;
}

    std::vector<ADB> BlackoilPropsAd::capPress(const ADB& sw,
                                               const ADB& so,
                                               const ADB& sg,
                                               const Cells& cells) const

    {
        const int numCells = cells.size();
        const int numActivePhases = numPhases();
        const int numBlocks = so.numBlocks();

        Block activeSat(numCells, numActivePhases);
        if (pu_.phase_used[Water]) {
            assert(sw.value().size() == numCells);
            activeSat.col(pu_.phase_pos[Water]) = sw.value();
        }
        if (pu_.phase_used[Oil]) {
            assert(so.value().size() == numCells);
            activeSat.col(pu_.phase_pos[Oil]) = so.value();
        } else {
            OPM_THROW(std::runtime_error, "BlackoilPropsAdFromDeck::relperm() assumes oil phase is active.");
        }
        if (pu_.phase_used[Gas]) {
            assert(sg.value().size() == numCells);
            activeSat.col(pu_.phase_pos[Gas]) = sg.value();
        }

        Block pc(numCells, numActivePhases);
        Block dpc(numCells, numActivePhases*numActivePhases);
        props_.capPress(numCells, activeSat.data(), cells.data(), pc.data(), dpc.data());

        std::vector<ADB> adbCapPressures;
        adbCapPressures.reserve(3);
        const ADB* s[3] = { &sw, &so, &sg };
        for (int phase1 = 0; phase1 < 3; ++phase1) {
            if (pu_.phase_used[phase1]) {
                const int phase1_pos = pu_.phase_pos[phase1];
                std::vector<ADB::M> jacs(numBlocks);
                for (int block = 0; block < numBlocks; ++block) {
                    jacs[block] = ADB::M(numCells, s[phase1]->derivative()[block].cols());
                }
                for (int phase2 = 0; phase2 < 3; ++phase2) {
                    if (!pu_.phase_used[phase2])
                        continue;
                    const int phase2_pos = pu_.phase_pos[phase2];
                    // Assemble dpc1/ds2.
                    const int column = phase1_pos + numActivePhases*phase2_pos; // Recall: Fortran ordering from props_.relperm()
                    ADB::M dpc1_ds2_diag = spdiag(dpc.col(column));
                    for (int block = 0; block < numBlocks; ++block) {
                        jacs[block] += dpc1_ds2_diag * s[phase2]->derivative()[block];
                    }
                }
                adbCapPressures.emplace_back(ADB::function(pc.col(phase1_pos), jacs));
            } else {
                adbCapPressures.emplace_back(ADB::null());
            }
        }
        return adbCapPressures;
    }

ConcreteProblem::ConcreteProblem(ConcreteProblemPtr parentProblemPtr, AbstractSolutionStepPtr solutionStepPtr) : AbstractProblem(parentProblemPtr->m_numberOfRow, parentProblemPtr->m_numberOfColumn) {
    m_freeCellPtrs = parentProblemPtr->m_freeCellPtrs;
    Cell cell1UsedByStep, cell2UsedByStep;
    solutionStepPtr->getCells(cell1UsedByStep, cell2UsedByStep);
    Cells cellsUsedByStep;
    cellsUsedByStep.push_back(cell1UsedByStep);
    cellsUsedByStep.push_back(cell2UsedByStep);
    regiterBlockedCells(cellsUsedByStep);
}

 auto valid_neighbours_of(Cell c) const
 {
     Cells neighbours;
     for (auto direction = '1'; direction != '9'; ++direction)
     {
         auto n = DIRECTIONS.at(direction)(c);
         if (validate(n))
             neighbours.insert(n);
     }
     return neighbours;
 }