C++ CI类代码示例

void compare_pressure(const GI& g, const std::vector<double>& p)
{
    typedef typename GI::CellIterator CI;
    int count = 0;
    double l1err = 0.0;
    double l2err = 0.0;
    double linferr = 0.0;
    double totv = 0.0;
    for (CI c = g.cellbegin(); c != g.cellend(); ++c, ++count) {
        Vec cen = c->centroid();
        double uval = u(cen);
        double diff = uval - p[count];
        double v = c->volume();
        l1err += std::fabs(diff*v);
        l2err += diff*diff*v;
        linferr = std::max(std::fabs(diff), linferr);
        totv += v;
        // std::cout << cen[0] << ' ' << uval << ' ' << p[count] << std::endl;
    }
    l2err = std::sqrt(l2err);
    std::cout << "\n\n"
              << "\n     L1 error density: " << l1err/totv
              << "\n     L2 error density: " << l2err/totv
              << "\n     Linf error:       " << linferr << "\n\n\n";
}

void assign_src(const GI& g, std::vector<double>& src)
{
    typedef typename GI::CellIterator CI;
    int count = 0;
    for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
        src[count++] = -Lu(c->centroid())*c->volume();
    }
}

void test_flowsolver(const GI& g, const RI& r)
{
    typedef typename GI::CellIterator                              CI;
    typedef typename CI::FaceIterator                              FI;
    typedef Opm::BasicBoundaryConditions<true, false>                  FBC;
    typedef Opm::IncompFlowSolverHybrid<GI, RI, FBC,
                                         Opm::MimeticIPAnisoRelpermEvaluator> FlowSolver;

    FlowSolver solver;

    typedef Opm::FlowBC BC;
    FBC flow_bc(7);
    //flow_bc.flowCond(1) = BC(BC::Dirichlet, 1.0*Opm::unit::barsa);
    //flow_bc.flowCond(2) = BC(BC::Dirichlet, 0.0*Opm::unit::barsa);
    flow_bc.flowCond(5) = BC(BC::Dirichlet, 100.0*Opm::unit::barsa);

    typename CI::Vector gravity;
    gravity[0] = gravity[1] = 0.0;
    gravity[2] = Opm::unit::gravity;

    solver.init(g, r, gravity, flow_bc);
    // solver.printStats(std::cout);

    //solver.assembleStatic(g, r);
    //solver.printIP(std::cout);

    std::vector<double> src(g.numberOfCells(), 0.0);
    std::vector<double> sat(g.numberOfCells(), 0.0);
#if 0
    if (g.numberOfCells() > 1) {
        src[0]     = 1.0;
        src.back() = -1.0;
    }
#endif

    solver.solve(r, sat, flow_bc, src);

#if 0
    solver.printSystem("system");
    typedef typename FlowSolver::SolutionType FlowSolution;
    FlowSolution soln = solver.getSolution();
    std::cout << "Cell Pressure:\n" << std::scientific << std::setprecision(15);
    for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
        std::cout << '\t' << soln.pressure(c) << '\n';
    }

    std::cout << "Cell (Out) Fluxes:\n";
    std::cout << "flux = [\n";
    for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
        for (FI f = c->facebegin(); f != c->faceend(); ++f) {
            std::cout << soln.outflux(f) << ' ';
        }
        std::cout << "\b\n";
    }
    std::cout << "]\n";
#endif
}

void test_evaluator(const Interface& g)
{
    typedef typename Interface::CellIterator CI;
    typedef typename CI       ::FaceIterator FI;
    typedef typename CI       ::Scalar       Scalar;

    typedef Dune::SharedFortranMatrix FMat;

    std::cout << "Called test_evaluator()" << std::endl;

    std::vector<int> numf; numf.reserve(g.numberOfCells());
    int max_nf = -1;
    for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
        numf.push_back(0);
        int& nf = numf.back();

        for (FI f = c->facebegin(); f != c->faceend(); ++f)
            ++nf;

        max_nf = std::max(max_nf, nf);
    }

    typedef int DummyClass;
    Dune::MimeticIPEvaluator<Interface, DummyClass> ip(max_nf);

    // Set dummy permeability K=diag(10,1,...,1,0.1).
    std::vector<Scalar> perm(dim * dim, Scalar(0.0));
    Dune::SharedCMatrix K(dim, dim, &perm[0]);
    for (int i = 0; i < dim; ++i)
        K(i,i) = 1.0;
    K(0    ,0    ) *= 10.0;
    K(dim-1,dim-1) /= 10.0;

    // Storage for inverse ip.
    std::vector<Scalar> ip_store(max_nf * max_nf, Scalar(0.0));

    // Loop grid whilst building (and outputing) the inverse IP matrix.
    int count = 0;
    for (CI c = g.cellbegin(); c != g.cellend(); ++c, ++count) {
        FMat Binv(numf[count], numf[count], &ip_store[0]);

        ip.evaluate(c, K, Binv);

        std::cout << count << " -> Binv = [\n" << Binv << "]\n";
    }
}

void assign_bc(const GI& g, BCS& bcs)
{
    typedef Dune::FlowBC BC;
    typedef typename GI::CellIterator CI;
    typedef typename CI::FaceIterator FI;
    int max_bid = 0;
    for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
        for (FI f = c->facebegin(); f != c->faceend(); ++f) {
            int bid = f->boundaryId();
            if (bid > max_bid) {
                max_bid = bid;
                bcs.resize(bid + 1);
            }
            bcs.flowCond(bid) = BC(BC::Dirichlet, u(f->centroid()));
        }
    }
}

 double setSourceBlock(const Vector& low, const Vector& high, double density, bool is_injection)
 {
     typedef typename GI::CellIterator CI;
     // Iterate over all cells, if the centroid of a cell is in the
     // domain given, set a source term. Accumulate total source terms.
     double total_rate = 0.0;
     for (CI c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c) {
         if (isInside(low, high, c->centroid())) {
             int index = c->index();
             double rate = density*c->volume();
             if (!is_injection) {
                 rate = -rate;
             }
             total_rate += rate;
             this->injection_rates_.addElement(rate, index);
             this->injection_rates_psolver_[index] = rate;
         }
     }
     return total_rate;
 }

    inline void setupRegionBasedConditions(const Opm::parameter::ParameterGroup& param,
                                           const GridInterface& g,
                                           BCs& bcs)
    {
        // Extract region and pressure value for Dirichlet bcs.
        typedef typename GridInterface::Vector Vector;
        Vector low;
        low[0] = param.getDefault("dir_block_low_x", 0.0);
        low[1] = param.getDefault("dir_block_low_y", 0.0);
        low[2] = param.getDefault("dir_block_low_z", 0.0);
        Vector high;
        high[0] = param.getDefault("dir_block_high_x", 1.0);
        high[1] = param.getDefault("dir_block_high_y", 1.0);
        high[2] = param.getDefault("dir_block_high_z", 1.0);
        double dir_block_pressure = param.get<double>("dir_block_pressure");

        // Set flow conditions for that region.
        // For this to work correctly, unique boundary ids should be used,
        // otherwise conditions may spread outside the given region, to all
        // faces with the same bid as faces inside the region.
        typedef typename GridInterface::CellIterator CI;
        typedef typename CI::FaceIterator FI;
        int max_bid = 0;
        std::vector<int> dir_bids;
        for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
            for (FI f = c->facebegin(); f != c->faceend(); ++f) {
                int bid = f->boundaryId();
                max_bid = std::max(bid, max_bid);
                if (bid != 0 && isInside(low, high, f->centroid())) {
                    dir_bids.push_back(bid);
                }
            }
        }
        bcs.resize(max_bid + 1);
        for (std::vector<int>::const_iterator it = dir_bids.begin(); it != dir_bids.end(); ++it) {
            bcs.flowCond(*it) = FlowBC(FlowBC::Dirichlet, dir_block_pressure);
        }

        // Transport BCs are defaulted.
    }

void CCSVWriter::WriteChanges (std::ostream& os, const CCachedLogInfo& cache)
{
	// header

	os << "ID,Revision,Change,PathID,CopyFromRev,CopyFromPathID\n";

	// content

	const CRevisionIndex& revisions = cache.GetRevisions();
	const CRevisionInfoContainer& logInfo = cache.GetLogInfo();

	// ids will be added on-the-fly

	size_t id = 0;
	for ( revision_t revision = revisions.GetFirstRevision()
		, last = revisions.GetLastRevision()
		; revision < last
		; ++revision)
	{
		index_t index = revisions[revision];
        if (index == NO_INDEX)
            continue;

		typedef CRevisionInfoContainer::CChangesIterator CI;

		if (  logInfo.GetPresenceFlags (index) 
			& CRevisionInfoContainer::HAS_CHANGEDPATHS)
		{
			// we actually have a valid (possibly empty) change list 
			// for this revision

			for ( CI iter = logInfo.GetChangesBegin (index)
				, end = logInfo.GetChangesEnd (index)
				; iter != end
				; ++iter)
			{
				static const char actions[9] = "AM R   D";
				char change = actions [iter->GetAction() /4 - 1];

				os << id++ << ',' 
				   << revision << ','
				   << change << ','
				   << iter->GetPathID() << ',';

				if (iter->HasFromPath())
				{
					os << iter->GetFromRevision() << ','
					   << iter->GetFromPathID() 
					   << "\n";
				}
				else
				{
					os << -1 << ','
					   << -1 
					   << "\n";
				}
			}
		}
	}
}

 void append(const CI &ci,
             const String &I, const Subspace<Spin> &S,
             std::vector<Excitation> &V) {
     // this particular traversal is used to generate semi-ordered list
     int idx = (int)ci.template strings<Spin>()[I];
     size_t i = I.size();
     while (i--) {
         if (!I[i]) continue; // empty orbital
         for (size_t j = 0; j < I.size(); ++j) {
             if (I[j] && (i != j)) continue; // not an empty orbital 
             String J = I.swap(i, j);
             if (!ci.test(J)) continue;
             int jdx = (int)ci.template strings<Spin>()[J];
             if (!S.test(jdx)) continue; // string not in S
             Excitation t;
             t.sgn = (float)sgn(I, i, j);
             t.integral = (int)index(i,j);
             t.I = idx;
             t.J = jdx;
             V.push_back(t);
         }
     }
 }

        void buildFaceIndices()
        {
#ifdef VERBOSE
            std::cout << "Building unique face indices... " << std::flush;
	    Opm::time::StopWatch clock;
            clock.start();
#endif
            typedef CellIterator CI;
            typedef typename CI::FaceIterator FI;

            // We build the actual cell to face mapping in two passes.
            // [code mostly lifted from IncompFlowSolverHybrid::enumerateGridDof(),
            //  but with a twist: This code builds a mapping from cells in index
            //  order to unique face numbers, while the mapping built in the
            //  enumerateGridDof() method was ordered by cell iterator order]

            // Allocate and reserve structures.
            const int nc = numberOfCells();
            std::vector<int> cell(nc, -1);
            std::vector<int> num_faces(nc); // In index order.
            std::vector<int> fpos;     fpos  .reserve(nc + 1);
            std::vector<int> num_cf;   num_cf.reserve(nc); // In iterator order.
            std::vector<int> faces ;

            // First pass: enumerate internal faces.
            int cellno = 0; fpos.push_back(0);
            int tot_ncf = 0, tot_ncf2 = 0, max_ncf = 0;
            for (CI c = cellbegin(); c != cellend(); ++c, ++cellno) {
                const int c0 = c->index();
                ASSERT((0 <= c0) && (c0 < nc) && (cell[c0] == -1));
                cell[c0] = cellno;
                num_cf.push_back(0);
                int& ncf = num_cf.back();
                for (FI f = c->facebegin(); f != c-> faceend(); ++f) {
                    if (!f->boundary()) {
                        const int c1 = f->neighbourCellIndex();
                        ASSERT((0 <= c1) && (c1 < nc) && (c1 != c0));

                        if (cell[c1] == -1) {
                            // Previously undiscovered internal face.
                            faces.push_back(c1);
                        }
                    }
                    ++ncf;
                }
                num_faces[c0] = ncf;
                fpos.push_back(int(faces.size()));
                max_ncf  = std::max(max_ncf, ncf);
                tot_ncf  += ncf;
                tot_ncf2 += ncf * ncf;
            }
            ASSERT(cellno == nc);

            // Build cumulative face sizes enabling direct insertion of
            // face indices into cfdata later.
            std::vector<int> cumul_num_faces(numberOfCells() + 1);
            cumul_num_faces[0] = 0;
            std::partial_sum(num_faces.begin(), num_faces.end(), cumul_num_faces.begin() + 1);

            // Avoid (most) allocation(s) inside 'c' loop.
            std::vector<int>    l2g;
            l2g.reserve(max_ncf);
            std::vector<double> cfdata(tot_ncf);
            int total_num_faces = int(faces.size());

            // Second pass: build cell-to-face mapping, including boundary.
            typedef std::vector<int>::iterator VII;
            for (CI c = cellbegin(); c != cellend(); ++c) {
                const int c0 = c->index();
                ASSERT ((0 <=      c0 ) && (     c0  < nc) &&
                        (0 <= cell[c0]) && (cell[c0] < nc));
                const int ncf = num_cf[cell[c0]];
                l2g.resize(ncf, 0);
                for (FI f = c->facebegin(); f != c->faceend(); ++f) {
                    if (f->boundary()) {
                        // External, not counted before.  Add new face...
                        l2g[f->localIndex()] = total_num_faces++;
                    } else {
                        // Internal face.  Need to determine during
                        // traversal of which cell we discovered this
                        // face first, and extract the face number
                        // from the 'faces' table range of that cell.

                        // Note: std::find() below is potentially
                        // *VERY* expensive (e.g., large number of
                        // seeks in moderately sized data in case of
                        // faulted cells).
                        const int c1 = f->neighbourCellIndex();
                        ASSERT ((0 <=      c1 ) && (     c1  < nc) &&
                                (0 <= cell[c1]) && (cell[c1] < nc));

                        int t = c0, seek = c1;
                        if (cell[seek] < cell[t])
                            std::swap(t, seek);
                        int s = fpos[cell[t]], e = fpos[cell[t] + 1];
                        VII p = std::find(faces.begin() + s, faces.begin() + e, seek);
                        ASSERT(p != faces.begin() + e);
                        l2g[f->localIndex()] = p - faces.begin();
                    }
                }
//.........这里部分代码省略.........

void AlgorithmMIG::compute_preliminary_blocks() throw (Exception) {
	CI* ci = NULL;

	long double* w_values = NULL;
	long double w_values_sum = 0.0;

	double lower_ci = 0.0;
	double upper_ci = 0.0;

	preliminary_block* new_strong_pairs = NULL;

	n_preliminary_blocks = 0u;

	ci = CIFactory::create(ci_method, likelihood_density);
	ci->set_dbview(db);

	w_values = (long double*)malloc(db->n_markers * sizeof(long double));
	if (w_values == NULL) {
		throw Exception(__FILE__, __LINE__, "Error in memory allocation.");
	}

	for (unsigned int i = 0u; i < db->n_markers; ++i) {
		w_values[i] = 0.0;
	}

	for (unsigned int i = 1u; i < db->n_markers; ++i) {
		w_values_sum = 0.0;
		for (long int j = i - 1u; j >= 0; --j) {
			ci->get_CI(i, j, &lower_ci, &upper_ci);
			if (!ISNAN(lower_ci) && !ISNAN(upper_ci)) {
				if (((auxiliary::fcmp(lower_ci, pos_strong_pair_cl, EPSILON) >= 0) && (auxiliary::fcmp(upper_ci, pos_strong_pair_cu, EPSILON) >= 0)) ||
						((auxiliary::fcmp(lower_ci, neg_strong_pair_cu, EPSILON) <= 0) && (auxiliary::fcmp(upper_ci, neg_strong_pair_cl, EPSILON) <= 0))) {
					w_values_sum += strong_pair_weight;
					w_values[j] += w_values_sum;
					if (auxiliary::fcmp(w_values[j], 0.0, EPSILON) >= 0) {
						if (n_preliminary_blocks >= preliminary_blocks_size) {
							preliminary_blocks_size += PRELIMINARY_BLOCKS_SIZE_INCREMENT;
							new_strong_pairs = (preliminary_block*)realloc(preliminary_blocks, preliminary_blocks_size * sizeof(preliminary_block));
							if (new_strong_pairs == NULL) {
								delete ci;
								ci = NULL;

								free(w_values);
								w_values = NULL;

								throw Exception(__FILE__, __LINE__, "Error in memory reallocation.");
							}
							preliminary_blocks = new_strong_pairs;
							new_strong_pairs = NULL;
						}

						preliminary_blocks[n_preliminary_blocks].start = j;
						preliminary_blocks[n_preliminary_blocks].end = i;
						preliminary_blocks[n_preliminary_blocks].length_bp = db->positions[i] - db->positions[j];

						++n_preliminary_blocks;
					}
				} else if ((auxiliary::fcmp(lower_ci, neg_recomb_pair_cu, EPSILON) >= 0) && (auxiliary::fcmp(upper_ci, pos_recomb_pair_cu, EPSILON) <= 0)) {
					w_values_sum -= recomb_pair_weight;
					w_values[j] += w_values_sum;
				} else {
					w_values[j] += w_values_sum;
				}
			} else {
				w_values[j] += w_values_sum;
			}
		}
	}

	delete ci;
	ci = NULL;

	free(w_values);
	w_values = NULL;
}

void test_flowsolver(const GI& g, const RI& r)
{
    typedef typename GI::CellIterator                   CI;
    typedef typename CI::FaceIterator                   FI;
    typedef Opm::BasicBoundaryConditions<true, false>  FBC;
    typedef Opm::IncompFlowSolverHybrid<GI, RI, FBC,
                                         Opm::MimeticIPEvaluator> FlowSolver;

    FlowSolver solver;

    typedef Opm::FlowBC BC;
    FBC flow_bc(7);

#if !USE_ALUGRID
    flow_bc.flowCond(5) = BC(BC::Dirichlet, 100.0*Opm::unit::barsa);
    flow_bc.flowCond(6) = BC(BC::Dirichlet, 0.0*Opm::unit::barsa);
#endif

    typename CI::Vector gravity(0.0);
    // gravity[2] = Dune::unit::gravity;

    solver.init(g, r, gravity, flow_bc);

    std::vector<double> src(g.numberOfCells(), 0.0);
    std::vector<double> sat(g.numberOfCells(), 0.0);
//     if (g.numberOfCells() > 1) {
//         src[0]     = 1.0;
//         src.back() = -1.0;
//     }

    solver.solve(r, sat, flow_bc, src, 5e-9, 3, 1);

#if 1
    typedef typename FlowSolver::SolutionType FlowSolution;
    FlowSolution soln = solver.getSolution();

    std::vector<typename GI::Vector> cell_velocity;
    estimateCellVelocity(cell_velocity, g, soln);
    // Dune's vtk writer wants multi-component data to be flattened.
    std::vector<double> cell_velocity_flat(&*cell_velocity.front().begin(),
                                           &*cell_velocity.back().end());
    std::vector<double> cell_pressure;
    getCellPressure(cell_pressure, g, soln);

    Dune::VTKWriter<typename GI::GridType::LeafGridView> vtkwriter(g.grid().leafView());
    vtkwriter.addCellData(cell_velocity_flat, "velocity", dim);
    vtkwriter.addCellData(cell_pressure, "pressure");
    vtkwriter.write("testsolution-" + boost::lexical_cast<std::string>(0),
                    Dune::VTKOptions::ascii);
#else    
    solver.printSystem("system");
    typedef typename FlowSolver::SolutionType FlowSolution;
    FlowSolution soln = solver.getSolution();

    std::cout << "Cell Pressure:\n" << std::scientific << std::setprecision(15);
    for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
        std::cout << '\t' << soln.pressure(c) << '\n';
    }

    std::cout << "Cell (Out) Fluxes:\n";
    std::cout << "flux = [\n";
    for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
        for (FI f = c->facebegin(); f != c->faceend(); ++f) {
            std::cout << soln.outflux(f) << ' ';
        }
        std::cout << "\b\n";
    }
    std::cout << "]\n";
#endif
}

	void update(const string& key,double val) {
		int idx = -1;
		bool found = da.exactMatchExists(key.c_str(),key.size(),&idx);
		if(found) {
			data[idx] = val;
			return;
		}
		insert(key,val);
	}

	void insert(const string& key,double val) {
		da.update(key.c_str(),key.size(),data.size());
		data.push_back(val);
	}

	bool exists(const string& key,double &old) {
		int idx = -1;
		bool found = da.exactMatchExists(key.c_str(),key.size(),&idx);
		if(found) old = data[idx];
		return found;
	}