C++ System类代码示例

void IBFEPostProcessor::interpolateVariables(const double data_time)
{
    Pointer<PatchHierarchy<NDIM> > hierarchy = d_fe_data_manager->getPatchHierarchy();
    const std::pair<int, int> patch_level_range = d_fe_data_manager->getPatchLevels();
    const int coarsest_ln = patch_level_range.first;
    const int finest_ln = patch_level_range.second - 1;

    const unsigned int num_eulerian_vars = d_scalar_interp_var_systems.size();

    // Set up Eulerian scratch space and fill ghost cell values.
    Pointer<RefineAlgorithm<NDIM> > refine_alg = new RefineAlgorithm<NDIM>();
    Pointer<RefineOperator<NDIM> > refine_op = NULL;
    std::set<int> scratch_idxs;
    for (unsigned int k = 0; k < num_eulerian_vars; ++k)
    {
        int& data_idx = d_scalar_interp_data_idxs[k];
        int& scratch_idx = d_scalar_interp_scratch_idxs[k];
        if (data_idx < 0 || scratch_idx < 0)
        {
            TBOX_ASSERT(data_idx < 0 || scratch_idx < 0);
            VariableDatabase<NDIM>* var_db = VariableDatabase<NDIM>::getDatabase();
            Pointer<hier::Variable<NDIM> > data_var = d_scalar_interp_vars[k];
            Pointer<VariableContext> data_ctx = d_scalar_interp_ctxs[k];
            data_idx = var_db->mapVariableAndContextToIndex(data_var, data_ctx);
            TBOX_ASSERT(data_idx >= 0);
            Pointer<VariableContext> scratch_ctx = var_db->getContext(d_name + "::SCRATCH");
            const FEDataManager::InterpSpec& interp_spec = d_scalar_interp_specs[k];
            const int ghost_width = LEInteractor::getMinimumGhostWidth(interp_spec.kernel_fcn);
            scratch_idx =
                var_db->registerVariableAndContext(data_var, scratch_ctx, ghost_width);
            scratch_idxs.insert(scratch_idx);
        }
        refine_alg->registerRefine(scratch_idx, data_idx, scratch_idx, refine_op);
    }
    CartExtrapPhysBdryOp refine_phys_bdry_op(scratch_idxs);
    for (int ln = coarsest_ln; ln <= finest_ln; ++ln)
    {
        Pointer<PatchLevel<NDIM> > level = hierarchy->getPatchLevel(ln);
        for (unsigned int k = 0; k < num_eulerian_vars; ++k)
        {
            const int scratch_idx = d_scalar_interp_scratch_idxs[k];
            if (!level->checkAllocated(scratch_idx))
                level->allocatePatchData(scratch_idx, data_time);
        }
        refine_alg->createSchedule(level, ln - 1, hierarchy, &refine_phys_bdry_op)
            ->fillData(data_time);
    }

    // Interpolate variables.
    NumericVector<double>* X_ghost_vec =
        d_fe_data_manager->buildGhostedCoordsVector(/*localize_data*/ true);
    for (unsigned int k = 0; k < num_eulerian_vars; ++k)
    {
        System* system = d_scalar_interp_var_systems[k];
        const std::string& system_name = system->name();
        const int scratch_idx = d_scalar_interp_scratch_idxs[k];
        d_fe_data_manager->interp(scratch_idx,
                                  *system->solution,
                                  *X_ghost_vec,
                                  system_name,
                                  d_scalar_interp_specs[k]);
    }

    // Deallocate Eulerian scratch space.
    for (int ln = coarsest_ln; ln <= finest_ln; ++ln)
    {
        Pointer<PatchLevel<NDIM> > level = hierarchy->getPatchLevel(ln);
        for (unsigned int k = 0; k < num_eulerian_vars; ++k)
        {
            const int scratch_idx = d_scalar_interp_scratch_idxs[k];
            level->deallocatePatchData(scratch_idx);
        }
    }
    return;
} // interpolateVariables

void testGBVI(GBVIForce::NonbondedMethod gbviMethod, CustomGBForce::NonbondedMethod customGbviMethod, std::string molecule) {

    const int numMolecules = 1;
    const double boxSize   = 10.0;
    ReferencePlatform platform;

    GBVIForce*        gbvi = new GBVIForce();
    std::vector<Vec3> positions;

    // select molecule

    if( molecule == "Monomer" ){
        buildMonomer( gbvi, positions );
    } else if( molecule == "Dimer" ){
        buildDimer( gbvi, positions );
    } else {
        buildEthane( gbvi, positions );
    }

    int numParticles = gbvi->getNumParticles();
    System standardSystem;
    System customGbviSystem;
    for (int i = 0; i < numParticles; i++) {
        standardSystem.addParticle(1.0);
        customGbviSystem.addParticle(1.0);
    }
    standardSystem.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0.0, 0.0), Vec3(0.0, boxSize, 0.0), Vec3(0.0, 0.0, boxSize));
    customGbviSystem.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0.0, 0.0), Vec3(0.0, boxSize, 0.0), Vec3(0.0, 0.0, boxSize));
    gbvi->setCutoffDistance(2.0);

    // create customGbviForce GBVI force

    CustomGBForce* customGbviForce  = createCustomGBVI( gbvi->getSolventDielectric(), gbvi->getSoluteDielectric() );
    customGbviForce->setCutoffDistance(2.0);

    // load parameters from gbvi to customGbviForce

    loadGbviParameters( gbvi, customGbviForce );

    OpenMM_SFMT::SFMT sfmt;
    init_gen_rand(0, sfmt);

    vector<Vec3> velocities(numParticles);
    for (int ii = 0; ii < numParticles; ii++) {
        velocities[ii] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
    }
    gbvi->setNonbondedMethod(gbviMethod);
    customGbviForce->setNonbondedMethod(customGbviMethod);

    standardSystem.addForce(gbvi);
    customGbviSystem.addForce(customGbviForce);

    VerletIntegrator integrator1(0.01);
    VerletIntegrator integrator2(0.01);

    Context context1(standardSystem, integrator1, platform);
    context1.setPositions(positions);
    context1.setVelocities(velocities);
    State state1 = context1.getState(State::Forces | State::Energy);

    Context context2(customGbviSystem, integrator2, platform);
    context2.setPositions(positions);
    context2.setVelocities(velocities);
    State state2 = context2.getState(State::Forces | State::Energy);

    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-4);

    for (int i = 0; i < numParticles; i++) {
        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-4);
    }
}

void DFFWidget::reloadHighLevelFile()
{
	setCurrentGeometry(NULL);

	if (frameModel)
		delete frameModel;
	if (geomModel)
		delete geomModel;
	if (mesh)
		delete mesh;

	DFFLoader dff;

	File file = dfile->getFile();

	try {
		mesh = dff.loadMesh(file);
		setEnabled(true);
	} catch (DFFException ex) {
		System::getInstance()->log(LogEntry::error(QString(tr("Error opening the DFF file: %1"))
				.arg(ex.getMessage().get()), &ex));
		mesh = NULL;
		setEnabled(false);
		//ui.hlWidget->setEnabled(false);
		//ui.tabWidget->setTabEnabled(ui.tabWidget->indexOf(ui.hlWidget), false);
	}

	// Load the frame tree
	frameModel = new DFFFrameItemModel(mesh);
	ui.frameTree->setModel(frameModel);

	connect(ui.frameTree->selectionModel(), SIGNAL(currentChanged(const QModelIndex&, const QModelIndex&)),
			this, SLOT(frameSelected(const QModelIndex&, const QModelIndex&)));

	System* sys = System::getInstance();

	DFFMesh::GeometryIterator git;

	if (mesh) {
		geomModel = new DFFGeometryItemModel(mesh);
		ui.geometryTree->setModel(geomModel);

		renderWidget->displayMesh(mesh);

		connect(ui.geometryTree->selectionModel(), SIGNAL(currentChanged(const QModelIndex&, const QModelIndex&)),
				this, SLOT(geometryTreeItemSelected(const QModelIndex&, const QModelIndex&)));

		connect(geomModel, SIGNAL(geometryDisplayStateChanged(DFFGeometry*, bool)), this,
				SLOT(geometryDisplayStateChanged(DFFGeometry*, bool)));
		connect(geomModel, SIGNAL(geometryPartDisplayStateChanged(DFFGeometryPart*, bool)), this,
				SLOT(geometryPartDisplayStateChanged(DFFGeometryPart*, bool)));

		// Search for the mesh texture
		char* meshName = new char[file.getPath().getFileName().length() + 1];
		strtolower(meshName, file.getPath().getFileName().get());
		meshName[strlen(meshName)-4] = '\0';

		SystemQuery query("FindMeshTextures");
		query["meshName"] = meshName;

		QList<SystemQueryResult> results = sys->sendSystemQuery(query);

		QStringList texes;

		for (unsigned int i = 0 ; i < results.size() ; i++) {
			SystemQueryResult res = results[i];
			QStringList subTexes = res["textures"].toStringList();

			texes.append(subTexes);
		}

		for (QStringList::iterator it = texes.begin() ; it != texes.end() ; it++) {
			QString tex = *it;
			ui.texSourceBox->addItem(tex, false);
		}

		if (texes.size() > 0)
			texSrcChanged(0);

		delete[] meshName;
	}
}

void testMembrane() {
    const int numMolecules = 70;
    const int numParticles = numMolecules*2;
    const double boxSize = 10.0;
    ReferencePlatform platform;

    // Create a system with an implicit membrane.

    System system;
    for (int i = 0; i < numParticles; i++) {
        system.addParticle(1.0);
    }
    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0.0, 0.0), Vec3(0.0, boxSize, 0.0), Vec3(0.0, 0.0, boxSize));
    CustomGBForce* custom = new CustomGBForce();
    custom->setCutoffDistance(2.0);
    custom->addPerParticleParameter("q");
    custom->addPerParticleParameter("radius");
    custom->addPerParticleParameter("scale");
    custom->addGlobalParameter("thickness", 3);
    custom->addGlobalParameter("solventDielectric", 78.3);
    custom->addGlobalParameter("soluteDielectric", 1);
    custom->addComputedValue("Imol", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(1/U^2-1/L^2)*(r-sr2*sr2/r)+0.5*log(L/U)/r+C);"
                             "U=r+sr2;"
                             "C=2*(1/or1-1/L)*step(sr2-r-or1);"
                             "L=max(or1, D);"
                             "D=abs(r-sr2);"
                             "sr2 = scale2*or2;"
                             "or1 = radius1-0.009; or2 = radius2-0.009", CustomGBForce::ParticlePairNoExclusions);
    custom->addComputedValue("Imem", "(1/radius+2*log(2)/thickness)/(1+exp(7.2*(abs(z)+radius-0.5*thickness)))", CustomGBForce::SingleParticle);
    custom->addComputedValue("B", "1/(1/or-tanh(1*psi-0.8*psi^2+4.85*psi^3)/radius);"
                             "psi=max(Imol,Imem)*or; or=radius-0.009", CustomGBForce::SingleParticle);
    custom->addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6-0.5*138.935456*(1/soluteDielectric-1/solventDielectric)*q^2/B", CustomGBForce::SingleParticle);
    custom->addEnergyTerm("-138.935456*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
                          "f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce::ParticlePairNoExclusions);
    vector<Vec3> positions(numParticles);
    vector<Vec3> velocities(numParticles);
    OpenMM_SFMT::SFMT sfmt;
    init_gen_rand(0, sfmt);
    vector<double> params(3);
    for (int i = 0; i < numMolecules; i++) {
        if (i < numMolecules/2) {
            params[0] = 1.0;
            params[1] = 0.2;
            params[2] = 0.5;
            custom->addParticle(params);
            params[0] = -1.0;
            params[1] = 0.1;
            custom->addParticle(params);
        }
        else {
            params[0] = 1.0;
            params[1] = 0.2;
            params[2] = 0.8;
            custom->addParticle(params);
            params[0] = -1.0;
            params[1] = 0.1;
            custom->addParticle(params);
        }
        positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
        positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
        velocities[2*i] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
        velocities[2*i+1] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
    }
    system.addForce(custom);
    VerletIntegrator integrator(0.01);
    Context context(system, integrator, platform);
    context.setPositions(positions);
    context.setVelocities(velocities);
    State state = context.getState(State::Forces | State::Energy);
    const vector<Vec3>& forces = state.getForces();

    // Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.

    double norm = 0.0;
    for (int i = 0; i < (int) forces.size(); ++i)
        norm += forces[i].dot(forces[i]);
    norm = std::sqrt(norm);
    const double stepSize = 1e-3;
    double step = 0.5*stepSize/norm;
    vector<Vec3> positions2(numParticles), positions3(numParticles);
    for (int i = 0; i < (int) positions.size(); ++i) {
        Vec3 p = positions[i];
        Vec3 f = forces[i];
        positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
        positions3[i] = Vec3(p[0]+f[0]*step, p[1]+f[1]*step, p[2]+f[2]*step);
    }
    context.setPositions(positions2);
    State state2 = context.getState(State::Energy);
    context.setPositions(positions3);
    State state3 = context.getState(State::Energy);
    ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/stepSize, 1e-3);
}

void testExclusions() {
    ReferencePlatform platform;
    for (int i = 3; i < 4; i++) {
        System system;
        system.addParticle(1.0);
        system.addParticle(1.0);
        VerletIntegrator integrator(0.01);
        CustomGBForce* force = new CustomGBForce();
        force->addComputedValue("a", "r", i < 2 ? CustomGBForce::ParticlePair : CustomGBForce::ParticlePairNoExclusions);
        force->addEnergyTerm("a", CustomGBForce::SingleParticle);
        force->addEnergyTerm("(1+a1+a2)*r", i%2 == 0 ? CustomGBForce::ParticlePair : CustomGBForce::ParticlePairNoExclusions);
        force->addParticle(vector<double>());
        force->addParticle(vector<double>());
        force->addExclusion(0, 1);
        system.addForce(force);
        Context context(system, integrator, platform);
        vector<Vec3> positions(2);
        positions[0] = Vec3(0, 0, 0);
        positions[1] = Vec3(1, 0, 0);
        context.setPositions(positions);
        State state = context.getState(State::Forces | State::Energy);
        const vector<Vec3>& forces = state.getForces();
        double f, energy;
        switch (i)
        {
            case 0: // e = 0
                f = 0;
                energy = 0;
                break;
            case 1: // e = r
                f = 1;
                energy = 1;
                break;
            case 2: // e = 2r
                f = 2;
                energy = 2;
                break;
            case 3: // e = 3r + 2r^2
                f = 7;
                energy = 5;
                break;
            default:
                ASSERT(false);
        }
        ASSERT_EQUAL_VEC(Vec3(f, 0, 0), forces[0], 1e-4);
        ASSERT_EQUAL_VEC(Vec3(-f, 0, 0), forces[1], 1e-4);
        ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), 1e-4);

        // Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.

        double norm = 0.0;
        for (int i = 0; i < (int) forces.size(); ++i)
            norm += forces[i].dot(forces[i]);
        norm = std::sqrt(norm);
        const double stepSize = 1e-3;
        double step = stepSize/norm;
        for (int i = 0; i < (int) positions.size(); ++i) {
            Vec3 p = positions[i];
            Vec3 f = forces[i];
            positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
        }
        context.setPositions(positions);
        State state2 = context.getState(State::Energy);
        ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state.getPotentialEnergy())/stepSize, 1e-3*abs(state.getPotentialEnergy()));
    }
}

void JumpErrorEstimator::estimate_error (const System& system,
                                         ErrorVector& error_per_cell,
                                         const NumericVector<Number>* solution_vector,
                                         bool estimate_parent_error)
{
  START_LOG("estimate_error()", "JumpErrorEstimator");
  /*

    Conventions for assigning the direction of the normal:

    - e & f are global element ids

    Case (1.) Elements are at the same level, e<f
    Compute the flux jump on the face and
    add it as a contribution to error_per_cell[e]
    and error_per_cell[f]

    ----------------------
    |           |          |
    |           |    f     |
    |           |          |
    |    e      |---> n    |
    |           |          |
    |           |          |
    ----------------------


    Case (2.) The neighbor is at a higher level.
    Compute the flux jump on e's face and
    add it as a contribution to error_per_cell[e]
    and error_per_cell[f]

    ----------------------
    |     |     |          |
    |     |  e  |---> n    |
    |     |     |          |
    |-----------|    f     |
    |     |     |          |
    |     |     |          |
    |     |     |          |
    ----------------------
  */

  // The current mesh
  const MeshBase& mesh = system.get_mesh();

  // The number of variables in the system
  const unsigned int n_vars = system.n_vars();

  // The DofMap for this system
  const DofMap& dof_map = system.get_dof_map();

  // Resize the error_per_cell vector to be
  // the number of elements, initialize it to 0.
  error_per_cell.resize (mesh.max_elem_id());
  std::fill (error_per_cell.begin(), error_per_cell.end(), 0.);

  // Declare a vector of floats which is as long as
  // error_per_cell above, and fill with zeros.  This vector will be
  // used to keep track of the number of edges (faces) on each active
  // element which are either:
  // 1) an internal edge
  // 2) an edge on a Neumann boundary for which a boundary condition
  //    function has been specified.
  // The error estimator can be scaled by the number of flux edges (faces)
  // which the element actually has to obtain a more uniform measure
  // of the error.  Use floats instead of ints since in case 2 (above)
  // f gets 1/2 of a flux face contribution from each of his
  // neighbors
  std::vector<float> n_flux_faces;
  if (scale_by_n_flux_faces)
    n_flux_faces.resize(error_per_cell.size(), 0);

  // Prepare current_local_solution to localize a non-standard
  // solution vector if necessary
  if (solution_vector && solution_vector != system.solution.get())
    {
      NumericVector<Number>* newsol =
        const_cast<NumericVector<Number>*>(solution_vector);
      System &sys = const_cast<System&>(system);
      newsol->swap(*sys.solution);
      sys.update();
    }

  fine_context.reset(new FEMContext(system));
  coarse_context.reset(new FEMContext(system));

  // Loop over all the variables we've been requested to find jumps in, to
  // pre-request
  for (var=0; var<n_vars; var++)
    {
      // Possibly skip this variable
      if (error_norm.weight(var) == 0.0) continue;

      // FIXME: Need to generalize this to vector-valued elements. [PB]
      FEBase* side_fe = NULL;

      const std::set<unsigned char>& elem_dims =
        fine_context->elem_dimensions();

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

void StatisticsSampler::sampleTemperature(System &system)
{
    m_temperature = 2.0/3.0 * m_kineticEnergy / system.atoms().size();
    system.setTemperature(m_temperature);
}

void RBEvaluation::read_in_vectors(System& sys,
                                   std::vector<NumericVector<Number>*>& vectors,
                                   const std::string& directory_name,
                                   const std::string& data_name,
                                   const bool read_binary_vectors)
{
  START_LOG("read_in_vectors()", "RBEvaluation");

  //libMesh::out << "Reading in the basis functions..." << std::endl;

  // Make sure processors are synced up before we begin
  this->comm().barrier();

  std::ostringstream file_name;
  const std::string basis_function_suffix = (read_binary_vectors ? ".xdr" : ".dat");
  struct stat stat_info;

  file_name << directory_name << "/" << data_name << "_header" << basis_function_suffix;
  Xdr header_data(file_name.str(),
                  read_binary_vectors ? DECODE : READ);

  // set the version number in header_data from io_version_string
  // (same code as in EquationSystemsIO::_read_impl)
  std::string io_version_string = get_io_version_string();
  std::string::size_type lm_pos = io_version_string.find("libMesh");
  std::istringstream iss(io_version_string.substr(lm_pos + 8));
  int ver_major = 0, ver_minor = 0, ver_patch = 0;
  char dot;
  iss >> ver_major >> dot >> ver_minor >> dot >> ver_patch;
  header_data.set_version(LIBMESH_VERSION_ID(ver_major, ver_minor, ver_patch));

  // We need to call sys.read_header (e.g. to set _written_var_indices properly),
  // but by setting the read_header argument to false, it doesn't reinitialize the system
  sys.read_header(header_data, io_version_string, /*read_header=*/false, /*read_additional_data=*/false);

  // Following EquationSystemsIO::read, we use a temporary numbering (node major)
  // before writing out the data
  MeshTools::Private::globally_renumber_nodes_and_elements(sys.get_mesh());


  const bool read_legacy_format = false;
  if (read_legacy_format)
    {
      // Use System::read_serialized_data to read in the basis functions
      // into this->solution and then swap with the appropriate
      // of basis function.
      for(unsigned int i=0; i<vectors.size(); i++)
	{
	  file_name.str(""); // reset the string
	  file_name << directory_name << "/" << data_name << i << basis_function_suffix;

	  // On processor zero check to be sure the file exists
	  if (this->processor_id() == 0)
	    {
	      int stat_result = stat(file_name.str().c_str(), &stat_info);

	      if (stat_result != 0)
		{
		  libMesh::out << "File does not exist: " << file_name.str() << std::endl;
		  libmesh_error();
		}
	    }

	  Xdr vector_data(file_name.str(),
		          read_binary_vectors ? DECODE : READ);

	  // The bf_data needs to know which version to read.
	  vector_data.set_version(LIBMESH_VERSION_ID(ver_major, ver_minor, ver_patch));

	  sys.read_serialized_data(vector_data, false);

	  vectors[i] = NumericVector<Number>::build(sys.comm()).release();
	  vectors[i]->init (sys.n_dofs(), sys.n_local_dofs(), false, libMeshEnums::PARALLEL);

	  // No need to copy, just swap
	  // *vectors[i] = *solution;
	  vectors[i]->swap(*sys.solution);
	}
    }

  //------------------------------------------------------
  // new implementation
  else
    {
      // Allocate storage for each vector
      for(unsigned int i=0; i<vectors.size(); i++)
	{
	  vectors[i] = NumericVector<Number>::build(sys.comm()).release();
	  vectors[i]->init (sys.n_dofs(), sys.n_local_dofs(), false, libMeshEnums::PARALLEL);
	}

      file_name.str("");
      file_name << directory_name << "/" << data_name << "_data" << basis_function_suffix;

      // On processor zero check to be sure the file exists
      if (this->processor_id() == 0)
	{
	  int stat_result = stat(file_name.str().c_str(), &stat_info);

	  if (stat_result != 0)
//.........这里部分代码省略.........

/*!
 *  \brief Initialize the operator
 *  \param ioSystem Current system
 */
void HPC::SendDemeToSupervisorOp::init(System& ioSystem)
{
    Beagle_StackTraceBeginM();
    mComm = castHandleT<HPC::MPICommunication>(ioSystem.getComponent("MPICommunication"));
    Beagle_StackTraceEndM("HPC::SendDemeToSupervisorOp::init(System&)");
}

int main(int argc, char* argv[])
{
	int program = 0;

	// Get input file from first command line argument
	if(argc < 2){ // No input file given
		std::cerr << "Usage: ./main [input_file]\n";
		program = -1;
	} else {
		std::string ifname = argv[1]; // Input filename
		std::string ofname = ifname; // Output file prefix
		std::size_t pos = ofname.find('.');
		if (pos != std::string::npos) { // Cut off extension
			ofname.erase(pos, ofname.length());
		}

		// Open the input file
		std::ifstream input(ifname);
		// Check it opened successfully
		if (!input.is_open()){
			std::cerr << "Failed to open input file.\n";
			program = -1;
		} else {
			
			// Make the system
			System sys = makeSystem(input);

			// Calculate the overlap integrals
			sys.calcOverlap();

			// Open main output file and print system details
			std::ofstream output(ofname + ".out");
			printSystem(sys, output, true);

			// Do all the optional commands
			int lastcmd = 0;
			int flag = 1;
			int orthog = 0;
			std::vector<int> currcmd;
			Eigen::MatrixXd f;
			while(flag > 0){
				currcmd = getNextCmd(input, lastcmd);
				switch(currcmd[0]){
				case 1: { // Print the overlap integrals
					std::ofstream intout(ofname + ".ints");
					printIntegrals(sys, intout);
					intout.close();
					break;
				}
				case 2: { // Print the sparse graph data
					std::ofstream sparseout(ofname + ".sparse");
					printSparseGraph(sys, sparseout, currcmd[1]);
					sparseout.close();
					break;
				}
				case 3: { // Canonical orthogonalisation
					f = orthogonalise(sys, currcmd[1], CANONICAL);
					orthog = 1;
					break;
				}
				case 4: { // Gram-Schmidt orthogonalisation
					f = orthogonalise(sys, currcmd[1], GRAM_SCHMIDT);
					orthog = 2;
					break;
				}
				case 5: { // Symmetric Lowdin orthogonalisation
					f = orthogonalise(sys, currcmd[1], SYM_LOWDIN); 
					orthog = 3;
					break;
				}
				case -1: { // Error
					output << "\nErroneous command given.\n";
					flag = 0;
					program = -1;
					break;
				}
				default: { // No more commands
					output << "\nProgram finished.\n";
					flag = 0;
				}
				}
				lastcmd++;
			}

			// Print orthogonalisation data if needed
			if (orthog > 0) {
				std::ofstream orthogout(ofname + ".orthog");
			    printOrthog(sys, orthogout, f, orthog);
				orthogout.close();
			}
			
			output.close();
		}
	}
	
	return program;
}

using namespace framework;

class Mode;
class View;
class Program;

Mode *		currentMode		= NULL;
View *		currentView		= NULL;
Mesh *		gridMesh		= NULL;
Mesh *		fullScreenQuad	= NULL;
Texture *	diffuseMap		= NULL,
                *	heightMap		= NULL,
                         *	normalMap		= NULL;
Program *	program			= NULL;
Query *		primitivesQuery	= NULL;
System		sys;
Input &		input			= *sys.getInput();
int			gridSize		= 127;
int			modeIndex		= 0;
int			currPart		= 0;
FPSCamera	cam;
bool		wireframe		= false;
bool		updateQuery		= true;
float		camSpeed		= 0.025f;

#include "program.h"
#include "modes.h"

int			viewIndex  = VIEW_COLOR;
Mode *		allModes[] = {	new Part1(), new Part2(), new Part3(),
              new Part4(), new Part5(), new Part6(), new Part7()

int main(int argc, char** argv)
{
	int round1, index;
	float round2;
	System *recSys;

	sys = new System* [10 * 11];

	for(round1 = 1; round1 < 11; round1++) {
		for(round2 = 1.00; round2 < 1.10; round2+=0.01) {
			index = ((round1-1) * 11) + ((round2 - 1.00) / 0.01);
			sys[index] = new System;

			#ifdef debug3
			printf("Round %d\n", round1+1);
			printf("File name: %s\n", argv[1]);
			#endif
			
			readInputFile(argv[1], sys[index]);
			setPrecedence(sys[index]);
			TAMwidthAssign(round1, round2, sys[index]);
			//sys[index].printPrecedence();
			sys[index]->setWaitExtList();
			sys[index]->setWaitBistList();
			sys[index]->TAMStat.initTAM(sys[index]->getSysTW());

			#ifdef debug1
			sys[index]->TAMStat.printTAM();
			printf("External List Size: %lu\n", sys[index]->ext_list.size());
			printf("Wiat External List Size: %lu\n", sys[index]->wait_ext_list.size());
			#endif

			sys[index]->fillTest(); //test this function
			//sys[index]->printResult(argv[1]);

			#ifdef debug2
			printf("Power Limit: %d\n", sys[index]->getSysPower());
			sys[index]->printExtList();
			sys[index]->printBistList();
			//sys[index]->TAMStat.printPowerStat();
			#endif 

			///*
			if(round1 == 1 && round2 == 1.00)
				recSys = sys[index];
			else
				recSys = getFasterSys(sys[index], recSys);
			//*/ 

			//delete sys;

			#ifdef debug3
			printf("Round1 %d & Round2 %f End!\n", round1, round2);
			#endif
		}
	}
	recSys->printResult(argv[1]);

	for(index = 0; index < 110; index++)
		delete sys[index];

	delete [] sys;

	return 0;
}

void ContinuousPanel::paint(const QRect& /*r*/, QPainter& p)
      {
      if (!_active) {
            _visible = false;
            return;
            }

      Measure* measure = _score->tick2measure(0);
      if (measure == 0){
            _visible = false;
            return;
            }

      if (measure->mmRest()) {
            measure = measure->mmRest();
            _mmRestCount = measure->mmRestCount();
            }

      System* system = measure->system();
            if (system == 0) {
                  _visible = false;
                  return;
                  }

      Segment* s = _score->tick2segment(0);
      double _spatium = _score->spatium();
      _x = 0;
      if (_width <= 0)
            _width  = s->x();

      //
      // Set panel height for whole system
      //
      _height = 6 * _spatium;
      _y = system->staffYpage(0) + system->page()->pos().y();
      double y2 = 0.0;
      for (int i = 0; i < _score->nstaves(); ++i) {
            SysStaff* ss = system->staff(i);
            if (!ss->show() || !_score->staff(i)->show())
                  continue;
            y2 = ss->y() + ss->bbox().height();
            }
      _height += y2 + 6*_spatium;
      _y -= 6 * _spatium;

      //
      // Check elements at current panel position
      //
      _offsetPanel = -(_sv->xoffset()) / _sv->mag();
      _rect = QRect(_offsetPanel + _width, _y, 1, _height);
      //qDebug() << "width=" << _width << "_y="<< _y << "_offsetPanel=" << _offsetPanel << "_sv->xoffset()" << _sv->xoffset() << "_sv->mag()" << _sv->mag() <<"_spatium" << _spatium << "s->canvasPos().x()" << s->canvasPos().x() << "s->x()" << s->x();
      Page* page = _score->pages().front();
      QList<Element*> elementsCurrent = page->items(_rect);
      if (elementsCurrent.empty()) {
            _visible = false;
            return;
            }
      qStableSort(elementsCurrent.begin(), elementsCurrent.end(), elementLessThan);

      _currentMeasure = nullptr;
      for (const Element* e : elementsCurrent) {
            e->itemDiscovered = 0;
            if (!e->visible()) {
                  if (_score->printing() || !_score->showInvisible())
                        continue;
                  }

            if (e->type() == Element::Type::MEASURE) {
                  _currentMeasure = static_cast<const Measure*>(e);
                  _currentTimeSig = _currentMeasure->timesig();
                  _currentMeasureTick = _currentMeasure->tick();
                  _currentMeasureNo = _currentMeasure->no();

                  // Find number of multi measure rests to display in the panel
                  if (_currentMeasure->isMMRest())
                        _mmRestCount = _currentMeasure->mmRestCount();
                  else if (_currentMeasure->mmRest())
                        _mmRestCount = _currentMeasure->mmRest()->mmRestCount();
                  else
                        _mmRestCount = 0;
                  _xPosMeasure = e->canvasX();
                  _measureWidth = e->width();
                  break;
                  }
            }
      if (_currentMeasure == nullptr)
            return;

      findElementWidths(elementsCurrent);

      // Don't show panel if staff names are visible
      if (_sv->xoffset() / _sv->mag() + _xPosMeasure > 0) {
            _visible = false;
            return;
            }
      //qDebug() << "_sv->xoffset()=" <<_sv->xoffset() << " _sv->mag()="<< _sv->mag() <<" s->x=" << s->x() << " width=" << _width << " currentMeasue=" << _currentMeasure->x() << " _xPosMeasure=" << _xPosMeasure;

      draw(p, elementsCurrent);
      _visible = true;
      }

void AdjointResidualErrorEstimator::estimate_error (const System& _system,
						    ErrorVector& error_per_cell,
						    const NumericVector<Number>* solution_vector,
					            bool estimate_parent_error)
{
  START_LOG("estimate_error()", "AdjointResidualErrorEstimator");

  // The current mesh
  const MeshBase& mesh = _system.get_mesh();

  // Resize the error_per_cell vector to be
  // the number of elements, initialize it to 0.
  error_per_cell.resize (mesh.max_elem_id());
  std::fill (error_per_cell.begin(), error_per_cell.end(), 0.);

  // Get the number of variables in the system
  unsigned int n_vars = _system.n_vars();

  // We need to make a map of the pointer to the solution vector
  std::map<const System*, const NumericVector<Number>*>solutionvecs;
  solutionvecs[&_system] = _system.solution.get();

  // Solve the dual problem if we have to
  if (!adjoint_already_solved)
    {
      // FIXME - we'll need to change a lot of APIs to make this trick
      // work with a const System...
      System&  system = const_cast<System&>(_system);
      system.adjoint_solve(_qoi_set);
    }

  // Flag to check whether we have not been asked to weight the variable error contributions in any specific manner
  bool error_norm_is_identity = error_norm.is_identity();

  // Create an ErrorMap/ErrorVector to store the primal, dual and total_dual variable errors
  ErrorMap primal_errors_per_cell;
  ErrorMap dual_errors_per_cell;
  ErrorMap total_dual_errors_per_cell;
  // Allocate ErrorVectors to this map if we're going to use it
  if (!error_norm_is_identity)
    for(unsigned int v = 0; v < n_vars; v++)
      {
        primal_errors_per_cell[std::make_pair(&_system, v)] = new ErrorVector;
        dual_errors_per_cell[std::make_pair(&_system, v)] = new ErrorVector;
        total_dual_errors_per_cell[std::make_pair(&_system, v)] = new ErrorVector;
      }
  ErrorVector primal_error_per_cell;
  ErrorVector dual_error_per_cell;
  ErrorVector total_dual_error_per_cell;

  // Have we been asked to weight the variable error contributions in any specific manner
  if(!error_norm_is_identity) // If we do
    {
      // Estimate the primal problem error for each variable
      _primal_error_estimator->estimate_errors
	(_system.get_equation_systems(), primal_errors_per_cell, &solutionvecs, estimate_parent_error);
    }
  else // If not
    {
      // Just get the combined error estimate
      _primal_error_estimator->estimate_error
	(_system, primal_error_per_cell, solution_vector, estimate_parent_error);
    }

  // Sum and weight the dual error estimate based on our QoISet
  for (unsigned int i = 0; i != _system.qoi.size(); ++i)
    {
      if (_qoi_set.has_index(i))
	{
	  // Get the weight for the current QoI
	  Real error_weight = _qoi_set.weight(i);

	   // We need to make a map of the pointer to the adjoint solution vector
	  std::map<const System*, const NumericVector<Number>*>adjointsolutionvecs;
	  adjointsolutionvecs[&_system] = &_system.get_adjoint_solution(i);

	  // Have we been asked to weight the variable error contributions in any specific manner
	  if(!error_norm_is_identity) // If we have
	    {
	      _dual_error_estimator->estimate_errors
		(_system.get_equation_systems(), dual_errors_per_cell, &adjointsolutionvecs,
		 estimate_parent_error);
	    }
 	  else // If not
	    {
	    // Just get the combined error estimate
	      _dual_error_estimator->estimate_error
		(_system, dual_error_per_cell, &(_system.get_adjoint_solution(i)), estimate_parent_error);
	    }

	  unsigned int error_size;

	  // Get the size of the first ErrorMap vector; this will give us the number of elements
	  if(!error_norm_is_identity) // If in non default weights case
            {
              error_size = dual_errors_per_cell[std::make_pair(&_system, 0)]->size();
            }
	  else // If in the standard default weights case
	    {
	      error_size = dual_error_per_cell.size();
//.........这里部分代码省略.........

CHECK(~ConformationSet())
	delete cs;
RESULT

CHECK(readDCDFile())
		ConformationSet cs;
		cs.readDCDFile(BALL_TEST_DATA_PATH(ConformationSet_test.dcd));
		cs.resetScoring();
		TEST_EQUAL(cs.size(), 10)
RESULT

CHECK(writeDCDFile(const String& filename, const Size num = 0))
		ConformationSet cs;

		PDBFile pdb(BALL_TEST_DATA_PATH(ConformationSet_test.pdb));
		System sys;
		pdb.read(sys);
		cs.setup(sys);

		cs.readDCDFile(BALL_TEST_DATA_PATH(ConformationSet_test.dcd));
		cs.resetScoring();
		String tmp_filename;
		NEW_TMP_FILE(tmp_filename)

		cs.writeDCDFile(tmp_filename);
		TEST_FILE(tmp_filename.c_str(), BALL_TEST_DATA_PATH(ConformationSet_test.dcd))
RESULT

CHECK(setup())
	PDBFile pdb(BALL_TEST_DATA_PATH(ConformationSet_test.pdb));
	System sys;