C++ ParameterList::validateParameters方法代码示例

void
ScatterScalarResponseBase<EvalT, Traits>::
setup(const Teuchos::ParameterList& p, const Teuchos::RCP<Albany::Layouts>& dl)
{
  bool stand_alone = p.get<bool>("Stand-alone Evaluator");

  // Setup fields we require
  PHX::Tag<ScalarT> global_response_tag =
    p.get<PHX::Tag<ScalarT> >("Global Response Field Tag");
  global_response = PHX::MDField<ScalarT>(global_response_tag);
  if (stand_alone)
    this->addDependentField(global_response);
  else
    this->addEvaluatedField(global_response);

  // Setup field we evaluate
  std::string fieldName = global_response_tag.name() + " Scatter Response";
  scatter_operation = Teuchos::rcp(new PHX::Tag<ScalarT>(fieldName, dl->dummy));
  this->addEvaluatedField(*scatter_operation);

  //! get and validate parameter list
  Teuchos::ParameterList* plist = 
    p.get<Teuchos::ParameterList*>("Parameter List");
  if (stand_alone) {
    Teuchos::RCP<const Teuchos::ParameterList> reflist = 
      this->getValidResponseParameters();
    plist->validateParameters(*reflist,0);
  }

  if (stand_alone)
    this->setName(fieldName+" Scatter Response");
}

TEUCHOS_UNIT_TEST(Ifpack2Parameters, Test0)
{
//we are now in a class method declared by the above macro, and
//that method has these input arguments:
//Teuchos::FancyOStream& out, bool& success

  Teuchos::ParameterList params;

  params.set("fact: iluk level-of-fill", (int) 2);

  Teuchos::ParameterList validparams;

  TEST_NOTHROW(Ifpack2::getValidParameters(validparams));

  params.validateParameters(validparams);

  int level_of_fill = 0;

  //call getParameter with a wrong name:
  Ifpack2::getParameter(params, "level-of-fill", level_of_fill);
  TEST_EQUALITY(level_of_fill, 0)

  //call getParameter with a valid name:
  Ifpack2::getParameter(params, "fact: iluk level-of-fill", level_of_fill);
  TEST_EQUALITY(level_of_fill, 2)
}

Aeras::ShallowWaterResponseL2Norm<EvalT, Traits>::
ShallowWaterResponseL2Norm(Teuchos::ParameterList& p,
		      const Teuchos::RCP<Albany::Layouts>& dl) :
  weighted_measure("Weights", dl->qp_scalar),
  flow_state_field("Flow State", dl->node_vector) 
{
  Teuchos::RCP<Teuchos::FancyOStream> out(Teuchos::VerboseObjectBase::getDefaultOStream());
  // get and validate Response parameter list
  Teuchos::ParameterList* plist = 
    p.get<Teuchos::ParameterList*>("Parameter List");
  Teuchos::RCP<const Teuchos::ParameterList> reflist = 
    this->getValidResponseParameters();
  plist->validateParameters(*reflist,0);

  std::string fieldName = "Flow Field"; //field to integral is the flow field 

  // coordinate dimensions
  std::vector<PHX::DataLayout::size_type> coord_dims;
  dl->qp_vector->dimensions(coord_dims);
  numQPs = coord_dims[1]; //# quad points
  numDims = coord_dims[2]; //# spatial dimensions
  std::vector<PHX::DataLayout::size_type> dims;
  flow_state_field.fieldTag().dataLayout().dimensions(dims);
  nPrimaryDOFs = 3; //we have 3 primary dofs for ShallowWater: h, u, v 
  numNodes =  dims[1]; //# nodes per element

 
  // User-specified parameters
  //None right now  

  // add dependent fields
  this->addDependentField(flow_state_field);
  this->addDependentField(weighted_measure);
  this->setName(fieldName+" Aeras Shallow Water L2 Norm"+PHX::typeAsString<EvalT>());
  
  using PHX::MDALayout;

  // Setup scatter evaluator
  p.set("Stand-alone Evaluator", false);
  std::string local_response_name = fieldName + " Local Response Aeras Shallow Water L2 Norm";
  std::string global_response_name = fieldName + " Global Response Aeras Shallow Water L2 Norm";
  int worksetSize = dl->qp_scalar->dimension(0);
  //There are four components of the response returned by this function: 
  //1.) |h|
  //2.) |u|
  //3.) |v|
  //4.) |solution vector|
  responseSize = 4; 
  Teuchos::RCP<PHX::DataLayout> local_response_layout = Teuchos::rcp(new MDALayout<Cell, Dim>(worksetSize, responseSize));
  Teuchos::RCP<PHX::DataLayout> global_response_layout = Teuchos::rcp(new MDALayout<Dim>(responseSize));
  PHX::Tag<ScalarT> local_response_tag(local_response_name, local_response_layout);
  PHX::Tag<ScalarT> global_response_tag(global_response_name, global_response_layout);
  p.set("Local Response Field Tag", local_response_tag);
  p.set("Global Response Field Tag", global_response_tag);
  PHAL::SeparableScatterScalarResponse<EvalT,Traits>::setup(p,dl);
}

int
Piro::PerformDakotaAnalysis(
    Thyra::ModelEvaluatorDefaultBase<double>& piroModel,
    Teuchos::ParameterList& dakotaParams,
    RCP< Thyra::VectorBase<double> >& p)
{
#ifdef HAVE_PIRO_TRIKOTA
  dakotaParams.validateParameters(*Piro::getValidPiroAnalysisDakotaParameters(),0);
  using std::string;

  string dakotaIn  = dakotaParams.get("Input File","dakota.in");
  string dakotaOut = dakotaParams.get("Output File","dakota.out");
  string dakotaErr = dakotaParams.get("Error File","dakota.err");
  string dakotaRes = dakotaParams.get("Restart File","dakota_restart.out");
  string dakotaRestartIn;
  if (dakotaParams.isParameter("Restart File To Read"))
    dakotaRestartIn = dakotaParams.get<string>("Restart File To Read");

  int dakotaRestartEvals= dakotaParams.get("Restart Evals To Read", 0);

  int p_index = dakotaParams.get("Parameter Vector Index", 0);
  int g_index = dakotaParams.get("Response Vector Index", 0);

  TriKota::Driver dakota(dakotaIn, dakotaOut, dakotaErr, dakotaRes,
                         dakotaRestartIn, dakotaRestartEvals);

  RCP<TriKota::ThyraDirectApplicInterface> trikota_interface =
    rcp(new TriKota::ThyraDirectApplicInterface
         (dakota.getProblemDescDB(), rcp(&piroModel,false), p_index, g_index),
	false);

  dakota.run(trikota_interface.get());

  Dakota::RealVector finalValues;
  if (dakota.rankZero())
    finalValues = dakota.getFinalSolution().all_continuous_variables();

  // Copy Dakota parameters into Thyra
  p = Thyra::createMember(piroModel.get_p_space(p_index));
  {
      Thyra::DetachedVectorView<double> global_p(p);
      for (int i = 0; i < finalValues.length(); ++i)
        global_p[i] = finalValues[i];
  }

  return 0;
#else
 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
 *out << "ERROR: Trilinos/Piro was not configured to include Dakota analysis."
      << "\nYou must enable TriKota." << endl;
 return 0;  // should not fail tests
#endif
}

QCAD::ResponseFieldValue<EvalT, Traits>::
ResponseFieldValue(Teuchos::ParameterList& p,
                   const Teuchos::RCP<Albany::Layouts>& dl) :
  coordVec("Coord Vec", dl->qp_vector),
  weights("Weights", dl->qp_scalar)
{
  // get and validate Response parameter list
  Teuchos::ParameterList* plist =
    p.get<Teuchos::ParameterList*>("Parameter List");
  Teuchos::RCP<const Teuchos::ParameterList> reflist =
    this->getValidResponseParameters();
  plist->validateParameters(*reflist,0);

  //! parameters passed down from problem
  Teuchos::RCP<Teuchos::ParameterList> paramsFromProblem =
    p.get< Teuchos::RCP<Teuchos::ParameterList> >("Parameters From Problem");

    // Material database (if given)
  if(paramsFromProblem != Teuchos::null)
    materialDB = paramsFromProblem->get< Teuchos::RCP<QCAD::MaterialDatabase> >("MaterialDB");
  else materialDB = Teuchos::null;

  // number of quad points per cell and dimension of space
  Teuchos::RCP<PHX::DataLayout> scalar_dl = dl->qp_scalar;
  Teuchos::RCP<PHX::DataLayout> vector_dl = dl->qp_vector;

  std::vector<PHX::DataLayout::size_type> dims;
  vector_dl->dimensions(dims);
  numQPs  = dims[1];
  numDims = dims[2];

  opRegion  = Teuchos::rcp( new QCAD::MeshRegion<EvalT, Traits>("Coord Vec","Weights",*plist,materialDB,dl) );

  // User-specified parameters
  operation    = plist->get<std::string>("Operation");

  bOpFieldIsVector = false;
  if(plist->isParameter("Operation Vector Field Name")) {
    opFieldName  = plist->get<std::string>("Operation Vector Field Name");
    bOpFieldIsVector = true;
  }
  else opFieldName  = plist->get<std::string>("Operation Field Name");

  bRetFieldIsVector = false;
  if(plist->isParameter("Return Vector Field Name")) {
    retFieldName  = plist->get<std::string>("Return Vector Field Name");
    bRetFieldIsVector = true;
  }
  else retFieldName = plist->get<std::string>("Return Field Name", opFieldName);
  bReturnOpField = (opFieldName == retFieldName);

  opX = plist->get<bool>("Operate on x-component", true) && (numDims > 0);
  opY = plist->get<bool>("Operate on y-component", true) && (numDims > 1);
  opZ = plist->get<bool>("Operate on z-component", true) && (numDims > 2);

  // setup operation field and return field (if it's a different field)
  if(bOpFieldIsVector) {
    PHX::MDField<ScalarT> f(opFieldName, vector_dl); opField = f; }
  else {
    PHX::MDField<ScalarT> f(opFieldName, scalar_dl); opField = f; }

  if(!bReturnOpField) {
    if(bRetFieldIsVector) {
      PHX::MDField<ScalarT> f(retFieldName, vector_dl); retField = f; }
    else {
      PHX::MDField<ScalarT> f(retFieldName, scalar_dl); retField = f; }
  }

  // add dependent fields
  this->addDependentField(opField);
  this->addDependentField(coordVec);
  this->addDependentField(weights);
  opRegion->addDependentFields(this);
  if(!bReturnOpField) this->addDependentField(retField); //when return field is *different* from op field

  // Set sentinal values for max/min problems
  initVals = Teuchos::Array<double>(5, 0.0);
  if( operation == "Maximize" ) initVals[1] = -1e200;
  else if( operation == "Minimize" ) initVals[1] = 1e100;
  else TEUCHOS_TEST_FOR_EXCEPTION (
    true, Teuchos::Exceptions::InvalidParameter, std::endl
    << "Error!  Invalid operation type " << operation << std::endl);

  this->setName(opFieldName+" Response Field Value"+PHX::TypeString<EvalT>::value);

  // Setup scatter evaluator
  std::string global_response_name =
    opFieldName + " Global Response Field Value";
  //int worksetSize = scalar_dl->dimension(0);
  int responseSize = 5;
  Teuchos::RCP<PHX::DataLayout> global_response_layout =
    Teuchos::rcp(new PHX::MDALayout<Dim>(responseSize));
  PHX::Tag<ScalarT> global_response_tag(global_response_name,
                                        global_response_layout);
  p.set("Stand-alone Evaluator", false);
  p.set("Global Response Field Tag", global_response_tag);
  this->setup(p,dl);

  // Specify which components of response (in this case 0th and 1st) to
  //  scatter derivatives for.
//.........这里部分代码省略.........

void InitialConditions(const Teuchos::RCP<Epetra_Vector>& soln,
                       const Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> > > >& wsElNodeEqID,
                       const Teuchos::ArrayRCP<std::string>& wsEBNames,
                       const Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > > coords,
                       const int neq, const int numDim,
                       Teuchos::ParameterList& icParams, const bool hasRestartSolution) {
  // Called twice, with x and xdot. Different param lists are sent in.
  icParams.validateParameters(*AAdapt::getValidInitialConditionParameters(wsEBNames), 0);

  // Default function is Constant, unless a Restart solution vector
  // was used, in which case the Init COnd defaults to Restart.
  std::string name;

  if(!hasRestartSolution) name = icParams.get("Function", "Constant");

  else                     name = icParams.get("Function", "Restart");

  if(name == "Restart") return;

  // Handle element block specific constant data

  if(name == "EBPerturb" || name == "EBPerturbGaussian" || name == "EBConstant") {

    bool perturb_values = false;

    Teuchos::Array<double> defaultData(neq);
    Teuchos::Array<double> perturb_mag;

    // Only perturb if the user has told us by how much to perturb
    if(name != "EBConstant" && icParams.isParameter("Perturb IC")) {

      perturb_values = true;

      perturb_mag = icParams.get("Perturb IC", defaultData);

    }

    /* The element block-based IC specification here is currently a hack. It assumes the initial value is constant
     * within each element across the element block (or optionally perturbed somewhat element by element). The
     * proper way to do this would be to project the element integration point values to the nodes using the basis
     * functions and a consistent mass matrix.
     *
     * The current implementation uses a single integration point per element - this integration point value for this
     * element within the element block is specified in the input file (and optionally perturbed). An approximation
     * of the load vector is obtained by accumulating the resulting (possibly perturbed) value into the nodes. Then,
     * a lumped version of the mass matrix is inverted and used to solve for the approximate nodal point initial
     * conditions.
     */

    // Use an Epetra_Vector to hold the lumped mass matrix (has entries only on the diagonal). Zero-ed out.

    Epetra_Vector lumpedMM(soln->Map(), true);

    // Make sure soln is zeroed - we are accumulating into it

    for(int i = 0; i < soln->MyLength(); i++)

      (*soln)[i] = 0;

    // Loop over all worksets, elements, all local nodes: compute soln as a function of coord and wsEBName


    Teuchos::RCP<AAdapt::AnalyticFunction> initFunc;

    for(int ws = 0; ws < wsElNodeEqID.size(); ws++) { // loop over worksets

      Teuchos::Array<double> data = icParams.get(wsEBNames[ws], defaultData);
      // Call factory method from library of initial condition functions

      if(perturb_values) {

        if(name == "EBPerturb")

          initFunc = Teuchos::rcp(new AAdapt::ConstantFunctionPerturbed(neq, numDim, ws, data, perturb_mag));

        else // name == EBGaussianPerturb

          initFunc = Teuchos::rcp(new
                                  AAdapt::ConstantFunctionGaussianPerturbed(neq, numDim, ws, data, perturb_mag));
      }

      else

        initFunc = Teuchos::rcp(new AAdapt::ConstantFunction(neq, numDim, data));

      std::vector<double> X(neq);
      std::vector<double> x(neq);

      for(int el = 0; el < wsElNodeEqID[ws].size(); el++) { // loop over elements in workset

        for(int i = 0; i < neq; i++)
          X[i] = 0;

        for(int ln = 0; ln < wsElNodeEqID[ws][el].size(); ln++) // loop over node local to the element
          for(int i = 0; i < neq; i++)
            X[i] += coords[ws][el][ln][i]; // nodal coords

        for(int i = 0; i < neq; i++)
          X[i] /= (double)neq;

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

Aeras::ShallowWaterResponseL2Error<EvalT, Traits>::
ShallowWaterResponseL2Error(Teuchos::ParameterList& p,
		      const Teuchos::RCP<Albany::Layouts>& dl) :
  sphere_coord("Lat-Long", dl->qp_gradient),
  weighted_measure("Weights", dl->qp_scalar),
  flow_state_field("Flow State", dl->node_vector), 
  BF("BF",dl->node_qp_scalar)
{
  Teuchos::RCP<Teuchos::FancyOStream> out(Teuchos::VerboseObjectBase::getDefaultOStream());
  // get and validate Response parameter list
  Teuchos::ParameterList* plist = 
    p.get<Teuchos::ParameterList*>("Parameter List");
  Teuchos::RCP<const Teuchos::ParameterList> reflist = 
    this->getValidResponseParameters();
  plist->validateParameters(*reflist,0);

  std::string fieldName = "Flow Field"; //field to integral is the flow field 

  // coordinate dimensions
  std::vector<PHX::DataLayout::size_type> coord_dims;
  dl->qp_vector->dimensions(coord_dims);
  numQPs = coord_dims[1]; //# quad points
  numDims = coord_dims[2]; //# spatial dimensions
  std::vector<PHX::DataLayout::size_type> dims;
  flow_state_field.fieldTag().dataLayout().dimensions(dims);
  vecDim = dims[2]; //# dofs per node
  numNodes =  dims[1]; //# nodes per element

 
  // User-specified parameters
  refSolName = plist->get<std::string>("Reference Solution Name"); //no reference solution by default.
  *out << "Reference Solution Name for Aeras::ShallowWaterResponseL2Error response: " << refSolName << std::endl; 
  inputData = plist->get<double>("Reference Solution Data", 0.0);
   
  if (refSolName == "Zero")
    ref_sol_name = ZERO;
  else if (refSolName == "TC2")
    ref_sol_name  = TC2;
  //Add other test case reference solutions here...
  else if (refSolName == "TC4"){
    ref_sol_name = TC4;
    
    myPi = Aeras::ShallowWaterConstants::self().pi;
    earthRadius = Aeras::ShallowWaterConstants::self().earthRadius;
    gravity = Aeras::ShallowWaterConstants::self().gravity;
    
    Omega = 2.0*myPi/(24.*3600.); //this should be sitting in SW Constants class
    
    rlon0 = 0.;
    rlat0 = myPi/4.;
    npwr = 14.;
    
    su0 = 20.;
    phi0 = 1.0e5;
    alfa = -0.03*(phi0/(2.*Omega*sin(myPi/4.)));
    sigma = (2.*earthRadius/1.0e6)*(2.*earthRadius/1.0e6);
    
  }
  else { 
    TEUCHOS_TEST_FOR_EXCEPTION(
      true, Teuchos::Exceptions::InvalidParameter,
      std::endl << "Error!  Unknown reference solution name " << ref_sol_name <<
      "!" << std::endl;);
  }

HydrostaticResponseL2Norm<EvalT, Traits>::
HydrostaticResponseL2Norm(Teuchos::ParameterList& p,
                     const Teuchos::RCP<Aeras::Layouts>& dl) :
  weighted_measure("Weights", dl->qp_scalar),
  velocity("Velx",  dl->qp_vector_level),
  temperature("Temperature",dl->qp_scalar_level),
  spressure("SPressure",dl->qp_scalar),
  numLevels(dl->node_scalar_level->dimension(2)), 
  out(Teuchos::VerboseObjectBase::getDefaultOStream())
{
  Teuchos::ParameterList* plist =
    p.get<Teuchos::ParameterList*>("Parameter List");
  Teuchos::RCP<const Teuchos::ParameterList> reflist =
    this->getValidResponseParameters();
  plist->validateParameters(*reflist,0);

  *out << "in Hydrostatic_Response_L2Norm! \n";

  // number of quad points per cell and dimension of space
  Teuchos::RCP<PHX::DataLayout> scalar_dl = dl->qp_scalar;
  Teuchos::RCP<PHX::DataLayout> vector_dl = dl->qp_vector;

  std::vector<PHX::DataLayout::size_type> dims;
  vector_dl->dimensions(dims);
  numQPs  = dims[1];
  numDims = dims[2];

  *out << "numQPs, numDims, numLevels: " << numQPs << ", " << numDims << ", " << numLevels << std::endl; 
  this->addDependentField(weighted_measure);
  this->addDependentField(velocity);
  this->addDependentField(temperature);
  this->addDependentField(spressure);

  this->setName("Aeras Hydrostatic Response L2 Norm");

  using PHX::MDALayout;

  // Setup scatter evaluator
  p.set("Stand-alone Evaluator", false);
  std::string local_response_name = "Local Response Aeras Hydrostatic Response L2 Norm";
  std::string global_response_name = "Global Response Aeras Hydrostatic Response L2 Norm";
  int worksetSize = scalar_dl->dimension(0);

  //FIXME: extend responseSize to have tracers 
  responseSize = 3*numLevels + 1; //there are 2 velocities and 1 temperature variable on each level
                                      //surface pressure is on 1st level only
                                      //the ordering is: Sp0, u0, v0, T0, u1, v1, T1, etc

  Teuchos::RCP<PHX::DataLayout> local_response_layout = Teuchos::rcp(
      new MDALayout<Cell,Dim>(worksetSize, responseSize));
  PHX::Tag<ScalarT> local_response_tag(local_response_name,
                                       local_response_layout);
  p.set("Local Response Field Tag", local_response_tag);

  Teuchos::RCP<PHX::DataLayout> global_response_layout = Teuchos::rcp(
      new MDALayout<Dim>(responseSize));
  PHX::Tag<ScalarT> global_response_tag(global_response_name,
                                        global_response_layout);
  p.set("Global Response Field Tag", global_response_tag);
  PHAL::SeparableScatterScalarResponse<EvalT,Traits>::setup(p,dl);

}

void
QCAD::PoissonProblem::constructDirichletEvaluators(
     const Albany::MeshSpecsStruct& meshSpecs)
{
   using Teuchos::RCP;
   using Teuchos::rcp;
   using Teuchos::ParameterList;
   using PHX::DataLayout;
   using PHX::MDALayout;
   using std::vector;
   using std::map;
   using std::string;

   using PHAL::DirichletFactoryTraits;
   using PHAL::AlbanyTraits;

   // Construct Dirichlet evaluators for all nodesets and names
   vector<string> dirichletNames(neq);
   dirichletNames[0] = "Phi";   
//   Albany::BCUtils<Albany::DirichletTraits> dirUtils;

   const std::vector<std::string>& nodeSetIDs = meshSpecs.nsNames;

   Teuchos::ParameterList DBCparams = params->sublist("Dirichlet BCs");
   DBCparams.validateParameters(*(Albany::DirichletTraits::getValidBCParameters(nodeSetIDs,dirichletNames)),0);

   map<string, RCP<ParameterList> > evaluators_to_build;
   RCP<DataLayout> dummy = rcp(new MDALayout<Dummy>(0));
   vector<string> dbcs;

   // Check for all possible standard BCs (every dof on every nodeset) to see which is set
   for (std::size_t i=0; i<nodeSetIDs.size(); i++) {
     for (std::size_t j=0; j<dirichletNames.size(); j++) {

       std::stringstream sstrm; sstrm << "DBC on NS " << nodeSetIDs[i] << " for DOF " << dirichletNames[j];
       std::string ss = sstrm.str();

       if (DBCparams.isParameter(ss)) {
         RCP<ParameterList> p = rcp(new ParameterList);
         int type = DirichletFactoryTraits<AlbanyTraits>::id_qcad_poisson_dirichlet;
         p->set<int>("Type", type);

         p->set< RCP<DataLayout> >("Data Layout", dummy);
         p->set< string >  ("Dirichlet Name", ss);
         p->set< RealType >("Dirichlet Value", DBCparams.get<double>(ss));
         p->set< string >  ("Node Set ID", nodeSetIDs[i]);
         p->set< int >     ("Number of Equations", dirichletNames.size());
         p->set< int >     ("Equation Offset", j);

         p->set<RCP<ParamLib> >("Parameter Library", paramLib);

         //! Additional parameters needed for Poisson Dirichlet BCs
         Teuchos::ParameterList& paramList = params->sublist("Poisson Source");
         p->set<Teuchos::ParameterList*>("Poisson Source Parameter List", &paramList);
         //p->set<string>("Temperature Name", "Temperature");  //to add if use shared param for DBC
         p->set<double>("Temperature", temperature);
         p->set< RCP<QCAD::MaterialDatabase> >("MaterialDB", materialDB);
	 p->set<double>("Energy unit in eV", energy_unit_in_eV);

         std::stringstream ess; ess << "Evaluator for " << ss;
         evaluators_to_build[ess.str()] = p;

         dbcs.push_back(ss);
       }
     }
   }

   //From here down, identical to Albany::AbstractProblem version of this function
   string allDBC="Evaluator for all Dirichlet BCs";
   {
      RCP<ParameterList> p = rcp(new ParameterList);
      int type = DirichletFactoryTraits<AlbanyTraits>::id_dirichlet_aggregator;
      p->set<int>("Type", type);

      p->set<vector<string>* >("DBC Names", &dbcs);
      p->set< RCP<DataLayout> >("Data Layout", dummy);
      p->set<string>("DBC Aggregator Name", allDBC);
      evaluators_to_build[allDBC] = p;
   }

   // Build Field Evaluators for each evaluation type
   PHX::EvaluatorFactory<AlbanyTraits,DirichletFactoryTraits<AlbanyTraits> > factory;
   RCP< vector< RCP<PHX::Evaluator_TemplateManager<AlbanyTraits> > > > evaluators;
   evaluators = factory.buildEvaluators(evaluators_to_build);

   // Create a DirichletFieldManager
   dfm = Teuchos::rcp(new PHX::FieldManager<AlbanyTraits>);

   // Register all Evaluators
   PHX::registerEvaluators(evaluators, *dfm);

   PHX::Tag<AlbanyTraits::Residual::ScalarT> res_tag0(allDBC, dummy);
   dfm->requireField<AlbanyTraits::Residual>(res_tag0);

   PHX::Tag<AlbanyTraits::Jacobian::ScalarT> jac_tag0(allDBC, dummy);
   dfm->requireField<AlbanyTraits::Jacobian>(jac_tag0);

   PHX::Tag<AlbanyTraits::Tangent::ScalarT> tan_tag0(allDBC, dummy);
   dfm->requireField<AlbanyTraits::Tangent>(tan_tag0);

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

PHAL::ResponseFieldIntegral<EvalT, Traits>::
ResponseFieldIntegral(Teuchos::ParameterList& p,
		      const Teuchos::RCP<Albany::Layouts>& dl) :
  coordVec("Coord Vec", dl->qp_gradient),
  weights("Weights", dl->qp_scalar)
{
  // get and validate Response parameter list
  Teuchos::ParameterList* plist = 
    p.get<Teuchos::ParameterList*>("Parameter List");
  Teuchos::RCP<const Teuchos::ParameterList> reflist = 
    this->getValidResponseParameters();
  plist->validateParameters(*reflist,0);

  // Get field type and corresponding layouts
  std::string field_name = plist->get<std::string>("Field Name");
  std::string fieldType = plist->get<std::string>("Field Type", "Scalar");
  if (plist->isType< Teuchos::Array<int> >("Field Components"))
    field_components = plist->get< Teuchos::Array<int> >("Field Components");
  Teuchos::RCP<PHX::DataLayout> field_layout;
  Teuchos::RCP<PHX::DataLayout> local_response_layout;
  Teuchos::RCP<PHX::DataLayout> global_response_layout;
  if (fieldType == "Scalar") {
    field_layout = dl->qp_scalar;
    local_response_layout = dl->cell_scalar;
    global_response_layout = dl->workset_scalar;
  }
  else if (fieldType == "Vector") {
    field_layout = dl->qp_vector;
    if (field_components.size() == 0) {
      local_response_layout = dl->cell_vector;
      global_response_layout = dl->workset_vector;
    }
    else {
      int worksetSize = dl->cell_scalar->dimension(0);
      local_response_layout = 
	Teuchos::rcp(new PHX::MDALayout<Cell,Dim>(worksetSize, 
						  field_components.size()));
      global_response_layout = 
	Teuchos::rcp(new PHX::MDALayout<Dim>(field_components.size()));
    }
  }
  else if (fieldType == "Tensor") {
    TEUCHOS_TEST_FOR_EXCEPTION(
      true, std::logic_error,
      "local_ and global_response must have rank 2. However, this code path "
      "makes them rank 3. Needs to be fixed.");
    field_layout = dl->qp_tensor;
    local_response_layout = dl->cell_tensor;
    global_response_layout = dl->workset_tensor;
  }
  else {
    TEUCHOS_TEST_FOR_EXCEPTION(
      true, 
      Teuchos::Exceptions::InvalidParameter,
      "Invalid field type " << fieldType << ".  Support values are " << 
      "Scalar, Vector, and Tensor." << std::endl);
  }
  field = PHX::MDField<ScalarT>(field_name, field_layout);
  field_layout->dimensions(field_dims);
  field_rank = field_layout->rank();
  if (field_components.size() == 0) {
    int num_components = field_dims[field_rank-1];
    field_components.resize(num_components);
    for (int i=0; i<num_components; i++)
      field_components[i] = i;
  }

  // coordinate dimensions
  std::vector<PHX::DataLayout::size_type> coord_dims;
  dl->qp_vector->dimensions(coord_dims);
  numQPs = coord_dims[1];
  numDims = coord_dims[2];
 
  // User-specified parameters
  std::string ebNameStr = plist->get<std::string>("Element Block Name","");
  if(ebNameStr.length() > 0) split(ebNameStr,',',ebNames);

  limitX = limitY = limitZ = false;
  if( plist->isParameter("x min") && plist->isParameter("x max") ) {
    limitX = true; TEUCHOS_TEST_FOR_EXCEPT(numDims <= 0);
    xmin = plist->get<double>("x min");
    xmax = plist->get<double>("x max");
  }
  if( plist->isParameter("y min") && plist->isParameter("y max") ) {
    limitY = true; TEUCHOS_TEST_FOR_EXCEPT(numDims <= 1);
    ymin = plist->get<double>("y min");
    ymax = plist->get<double>("y max");
  }
  if( plist->isParameter("z min") && plist->isParameter("z max") ) {
    limitZ = true; TEUCHOS_TEST_FOR_EXCEPT(numDims <= 2);
    zmin = plist->get<double>("z min");
    zmax = plist->get<double>("z max");
  }

  // length scaling
  double X0 = plist->get<double>("Length Scaling", 1.0);
  if (numDims == 1)
    scaling = X0; 
  else if (numDims == 2)
    scaling = X0*X0; 
//.........这里部分代码省略.........

int
Piro::PerformAnalysis(
    Thyra::ModelEvaluatorDefaultBase<double>& piroModel,
    Teuchos::ParameterList& analysisParams,
    RCP< Thyra::VectorBase<double> >& result)
{

  analysisParams.validateParameters(*Piro::getValidPiroAnalysisParameters(),0);

  int status;
  RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();

  string analysis = analysisParams.get<string>("Analysis Package");
  *out << "\n\nPiro::PerformAnalysis() requests: " << analysis << endl;

  if (analysis=="Solve") {
    *out << "Piro PerformAnalysis: Model Solve Being Performed " << endl;
    Piro::PerformSolveBase(piroModel, analysisParams.sublist("Solve"), result);
    status = 0; // Succeeds or throws
  }
#ifdef HAVE_PIRO_TRIKOTA
  else if (analysis=="Dakota") {
    *out << "Piro PerformAnalysis: Dakota Analysis Being Performed " << endl;

    status = Piro::PerformDakotaAnalysis(piroModel,
                         analysisParams.sublist("Dakota"), result);

  }
#endif
#ifdef HAVE_PIRO_MOOCHO
  else if (analysis == "MOOCHO") {
    *out << "Piro PerformAnalysis: MOOCHO Optimization Being Performed " << endl;
    status = Piro::PerformMoochoAnalysis(piroModel,
                          analysisParams.sublist("MOOCHO"), result);

  }
#endif
#ifdef HAVE_PIRO_OPTIPACK
  else if (analysis == "OptiPack") {
    *out << "Piro PerformAnalysis: Optipack Optimization Being Performed " << endl;
    status = Piro::PerformOptiPackAnalysis(piroModel,
                    analysisParams.sublist("OptiPack"),
                    analysisParams.sublist("GlobiPack"), result);

  }
#endif
#ifdef HAVE_PIRO_ROL
  else if (analysis == "ROL") {
    *out << "Piro PerformAnalysis: ROL Optimization Being Performed " << endl;
    status = Piro::PerformROLAnalysis(piroModel,
                          analysisParams.sublist("ROL"), result);

  }
#endif
  else {
    if (analysis == "Dakota" || analysis == "OptiPack" || analysis == "MOOCHO" || analysis == "ROL")
      *out << "ERROR: Trilinos/Piro was not configured to include \n "
           << "       analysis type: " << analysis << endl;
    else
      *out << "ERROR: Piro: Unknown analysis type: " << analysis << "\n"
           << "       Valid analysis types are: Solve, Dakota, MOOCHO, OptiPack, ROL\n" << endl;
    status = 0; // Should not fail tests
  }

  // Output status and paramters
  if (status==0)  *out << "\nPiro Analysis Finished successfully." << endl;
  else  *out << "\nPiro Analysis failed with status: " << status << endl;

  if ( analysisParams.get("Output Final Parameters", true) )
    if (result != Teuchos::null) {
       *out << "\tFinal parameters are: " << "\n\tp = ";
       *out << Teuchos::describe(*result, Teuchos::VERB_EXTREME ) << endl;
    }

  return status;
}

QCAD::ResponseSaddleValue<EvalT, Traits>::
ResponseSaddleValue(Teuchos::ParameterList& p,
		    const Teuchos::RCP<Albany::Layouts>& dl) :
  coordVec(p.get<std::string>("Coordinate Vector Name"), dl->qp_vector),
  coordVec_vertices(p.get<std::string>("Coordinate Vector Name"), dl->vertices_vector),
  weights(p.get<std::string>("Weights Name"), dl->qp_scalar)
{
  using Teuchos::RCP;
  
  //! get lattice temperature and materialDB from "Parameters From Problem"
  RCP<Teuchos::ParameterList> probList = 
    p.get< RCP<Teuchos::ParameterList> >("Parameters From Problem");
  lattTemp = probList->get<double>("Temperature");
  materialDB = probList->get< RCP<QCAD::MaterialDatabase> >("MaterialDB");
  
  //! get and validate Response parameter list
  Teuchos::ParameterList* plist = 
    p.get<Teuchos::ParameterList*>("Parameter List");

  Teuchos::RCP<const Teuchos::ParameterList> reflist = 
    this->getValidResponseParameters();
  plist->validateParameters(*reflist,0);

  //! get pointer to response function object
  svResponseFn = plist->get<Teuchos::RCP<QCAD::SaddleValueResponseFunction> >
    ("Response Function");

  //! number of quad points per cell and dimension of space
  std::vector<PHX::DataLayout::size_type> dims;
  dl->qp_vector->dimensions(dims);
  numQPs  = dims[1];
  numDims = dims[2];
  dl->vertices_vector->dimensions(dims);
  numVertices = dims[2];

  //! User-specified parameters
  fieldName  = plist->get<std::string>("Field Name");
  
  // limit to Potential only because other fields such as CB show large error 
  // and very jaggy profile, which may be due to averaging effect because Ec is not
  // well defined at the Si/SiO2 interface (discontinuous), while Potential is always continuous.
  if (fieldName != "Potential") 
     TEUCHOS_TEST_FOR_EXCEPTION (true, Teuchos::Exceptions::InvalidParameter, std::endl 
		      << "Error! Field Name must be Potential" << std::endl); 

  fieldGradientName  = plist->get<std::string>("Field Gradient Name");
  scaling = plist->get<double>("Field Scaling Factor",-1.0);
  gradScaling = plist->get<double>("Field Gradient Scaling Factor",-1.0);

  retFieldName = plist->get<std::string>("Return Field Name", fieldName);
  retScaling = plist->get<double>("Return Field Scaling Factor",1.0);
  bReturnSameField = (fieldName == retFieldName);
  //bLateralVolumes = true; // Future: make into a parameter

  //! Special case when return field name == "current": then just compute 
  //   as if returning the same field, and overwrite with current value at end
  if(retFieldName == "current")
    bReturnSameField = true;    

  //! setup operation field and its gradient, and the return field (if it's different)
  PHX::MDField<ScalarT> f(fieldName, dl->qp_scalar); field = f;
  PHX::MDField<ScalarT> fg(fieldGradientName, dl->qp_vector); fieldGradient = fg;

  if(!bReturnSameField) {
    PHX::MDField<ScalarT> fr(retFieldName, dl->qp_scalar); retField = fr; }


  //! add dependent fields
  this->addDependentField(field);
  this->addDependentField(fieldGradient);
  this->addDependentField(coordVec);
  this->addDependentField(coordVec_vertices);
  this->addDependentField(weights);
  if(!bReturnSameField) this->addDependentField(retField);

  std::string responseID = "QCAD Saddle Value";
  this->setName(responseID + PHX::TypeString<EvalT>::value);

  /*//! response evaluator must evaluate dummy operation
  Teuchos::RCP<PHX::DataLayout> dummy_dl =
    p.get< Teuchos::RCP<PHX::DataLayout> >("Dummy Data Layout");
  
  response_operation = Teuchos::rcp(new PHX::Tag<ScalarT>(responseID, dummy_dl));
  this->addEvaluatedField(*response_operation);*/

  // Setup scatter evaluator
  p.set("Stand-alone Evaluator", false);

  int responseSize = 5;
  int worksetSize = dl->qp_scalar->dimension(0);
  Teuchos::RCP<PHX::DataLayout> global_response_layout =
    Teuchos::rcp(new PHX::MDALayout<Dim>(responseSize));
  Teuchos::RCP<PHX::DataLayout> local_response_layout =
    Teuchos::rcp(new PHX::MDALayout<Cell,Dim>(worksetSize, responseSize));


  std::string local_response_name = 
    fieldName + " Local Response Saddle Value";
  std::string global_response_name = 
    fieldName + " Global Response Saddle Value";
//.........这里部分代码省略.........

Adapt::ElementSizeFieldBase<EvalT, Traits>::
ElementSizeFieldBase(Teuchos::ParameterList& p,
		  const Teuchos::RCP<Albany::Layouts>& dl) :
  coordVec(p.get<std::string>("Coordinate Vector Name"), dl->qp_vector),
  coordVec_vertices(p.get<std::string>("Coordinate Vector Name"), dl->vertices_vector),
  qp_weights("Weights", dl->qp_scalar)
{

  //! get and validate ElementSizeField parameter list
  Teuchos::ParameterList* plist =
    p.get<Teuchos::ParameterList*>("Parameter List");
  Teuchos::RCP<const Teuchos::ParameterList> reflist =
    this->getValidSizeFieldParameters();
  plist->validateParameters(*reflist,0);

  // Isotropic --> element size scalar corresponding to the nominal element radius
  // Anisotropic --> element size vector (x, y, z) with the width, length, and height of the element
  // Weighted versions (upcoming) --> scale the above sizes with a scalar or vector field

  className = plist->get<std::string>("Size Field Name", "Element_Size_Field");
  outputToExodus = plist->get<bool>("Output to File", true);
  outputCellAverage = plist->get<bool>("Generate Cell Average", true);
  outputQPData = plist->get<bool>("Generate QP Values", false);
  outputNodeData = plist->get<bool>("Generate Nodal Values", false);
  isAnisotropic = plist->get<bool>("Anisotropic Size Field", false);

  //! number of quad points per cell and dimension
  Teuchos::RCP<PHX::DataLayout> scalar_dl = dl->qp_scalar;
  Teuchos::RCP<PHX::DataLayout> vector_dl = dl->qp_vector;
  Teuchos::RCP<PHX::DataLayout> cell_dl = dl->cell_scalar;
  Teuchos::RCP<PHX::DataLayout> vert_vector_dl = dl->vertices_vector;
  numQPs = vector_dl->dimension(1);
  numDims = vector_dl->dimension(2);
  numVertices = vert_vector_dl->dimension(2);
  this->addDependentField(qp_weights);
  this->addDependentField(coordVec);
  this->addDependentField(coordVec_vertices);

  //! Register with state manager
  this->pStateMgr = p.get< Albany::StateManager* >("State Manager Ptr");

  if( outputCellAverage ) {
    if(isAnisotropic) //An-isotropic
      this->pStateMgr->registerStateVariable(className + "_Cell", dl->cell_vector, dl->dummy, "all", "scalar",
         0.0, false, outputToExodus);
    else
      this->pStateMgr->registerStateVariable(className + "_Cell", dl->cell_scalar, dl->dummy, "all", "scalar",
         0.0, false, outputToExodus);
  }

  if( outputQPData ) {
//    if(isAnisotropic) //An-isotropic
//    Always anisotropic?
      this->pStateMgr->registerStateVariable(className + "_QP", dl->qp_vector, dl->dummy, "all",
        "scalar", 0.0, false, outputToExodus);
//    else
//      this->pStateMgr->registerStateVariable(className + "_QP", dl->qp_scalar, dl->dummy, "all",
//        "scalar", 0.0, false, outputToExodus);
  }

  if( outputNodeData ) {
    // The weighted projected value

    // Note that all dl->node_node_* layouts are handled by the Adapt_NodalDataBlock class, inside
    // of the state manager, as they require interprocessor synchronization

    if(isAnisotropic){ //An-isotropic
      this->pStateMgr->registerNodalBlockStateVariable(className + "_Node", dl->node_node_vector, dl->dummy, "all",
         "scalar", 0.0, false, outputToExodus);

    }
    else {
      this->pStateMgr->registerNodalBlockStateVariable(className + "_Node", dl->node_node_scalar, dl->dummy, "all",
         "scalar", 0.0, false, outputToExodus);


    }

    // The value of the weights used in the projection
    // Initialize to zero - should give us nan's during the division step if something is wrong
    this->pStateMgr->registerNodalBlockStateVariable(className + "_NodeWgt", dl->node_node_scalar, dl->dummy, "all",
         "scalar", 0.0, false, outputToExodus);

  }

  // Create field tag
  size_field_tag =
    Teuchos::rcp(new PHX::Tag<ScalarT>(className, dl->dummy));

  this->addEvaluatedField(*size_field_tag);
}

TotalVolume<EvalT, Traits>::
TotalVolume(Teuchos::ParameterList& p,
                     const Teuchos::RCP<Aeras::Layouts>& dl) :
  coordVec("Coord Vec", dl->qp_vector),
  weighted_measure("Weights", dl->qp_scalar),
  density ("Density", dl->qp_scalar_level),
  velocity("Velx",  dl->qp_vector_level),
  temperature("Temperature",dl->qp_scalar_level),
  Cpstar("Cpstar",dl->qp_scalar_level),
  pie("Pi",  dl->qp_scalar_level),
   numLevels(dl->node_scalar_level->dimension(2))


{
  Teuchos::ParameterList* plist =
    p.get<Teuchos::ParameterList*>("Parameter List");
  Teuchos::RCP<const Teuchos::ParameterList> reflist =
    this->getValidResponseParameters();
  plist->validateParameters(*reflist,0);

  Phi0 = 0;
   std::cout << "Total_Volume: Phi0 = " << Phi0 << std::endl;

  // number of quad points per cell and dimension of space
  Teuchos::RCP<PHX::DataLayout> scalar_dl = dl->qp_scalar;
  Teuchos::RCP<PHX::DataLayout> vector_dl = dl->qp_vector;

  std::vector<PHX::DataLayout::size_type> dims;
  vector_dl->dimensions(dims);
  numQPs  = dims[1];
  numDims = dims[2];

//  this->addDependentField(coordVec);
  this->addDependentField(weighted_measure);
  this->addDependentField(density);
  this->addDependentField(velocity);
  this->addDependentField(temperature);
  this->addDependentField(Cpstar);
  this->addDependentField(pie);

  this->setName("Aeras Total Volume");

  using PHX::MDALayout;

  // Setup scatter evaluator
  p.set("Stand-alone Evaluator", false);
  std::string local_response_name = "Local Response Aeras Total Volume";
  std::string global_response_name = "Global Response Aeras Total Volume";
  int worksetSize = scalar_dl->dimension(0);

  int responseSize = 3;

  Teuchos::RCP<PHX::DataLayout> local_response_layout = Teuchos::rcp(
      new MDALayout<Cell,Dim>(worksetSize, responseSize));
  PHX::Tag<ScalarT> local_response_tag(local_response_name,
                                       local_response_layout);
  p.set("Local Response Field Tag", local_response_tag);

  Teuchos::RCP<PHX::DataLayout> global_response_layout = Teuchos::rcp(
      new MDALayout<Dim>(responseSize));
  PHX::Tag<ScalarT> global_response_tag(global_response_name,
                                        global_response_layout);
  p.set("Global Response Field Tag", global_response_tag);
  PHAL::SeparableScatterScalarResponse<EvalT,Traits>::setup(p,dl);

}

QCAD::ResponseFieldIntegral<EvalT, Traits>::
ResponseFieldIntegral(Teuchos::ParameterList& p,
		      const Teuchos::RCP<Albany::Layouts>& dl) :
  coordVec("Coord Vec", dl->qp_vector),
  weights("Weights", dl->qp_scalar)
{
  //! get and validate Response parameter list
  Teuchos::ParameterList* plist = 
    p.get<Teuchos::ParameterList*>("Parameter List");
  Teuchos::RCP<const Teuchos::ParameterList> reflist = 
    this->getValidResponseParameters();
  plist->validateParameters(*reflist,0);

  //! parameters passed down from problem
  Teuchos::RCP<Teuchos::ParameterList> paramsFromProblem = 
    p.get< Teuchos::RCP<Teuchos::ParameterList> >("Parameters From Problem");
  if(paramsFromProblem != Teuchos::null) {

    // Material database 
    materialDB = paramsFromProblem->get< Teuchos::RCP<QCAD::MaterialDatabase> >("MaterialDB");

    // Length unit in meters
    length_unit_in_m = paramsFromProblem->get<double>("Length unit in m");
  }
  else {
    materialDB = Teuchos::null;
    length_unit_in_m = 1.0e-6; //default length unit = microns (backward compat)
  }
       
  //! number of quad points per cell
  Teuchos::RCP<PHX::DataLayout> scalar_dl = dl->qp_scalar;
  numQPs = scalar_dl->dimension(1);
  
  //! obtain number of dimensions
  Teuchos::RCP<PHX::DataLayout> vector_dl = dl->qp_vector;
  std::vector<PHX::DataLayout::size_type> dims;
  vector_dl->dimensions(dims);
  numDims = dims[2];

  //! Initialize Region
  opRegion  = Teuchos::rcp( new QCAD::MeshRegion<EvalT, Traits>("Coord Vec","Weights",*plist,materialDB,dl) );

  //! User-specified parameters
  std::string fieldName;

  fieldName = plist->get<std::string>("Field Name","");
  if(fieldName.length() > 0) {
    fieldNames.push_back(fieldName);
    conjugateFieldFlag.push_back(plist->get<bool>("Conjugate Field",false));

    fieldName = plist->get<std::string>("Field Name Im","");
    fieldNames_Imag.push_back(fieldName);
    if(fieldName.length() > 0) fieldIsComplex.push_back(true);
    else fieldIsComplex.push_back(false);
  }

  for(int i=1; i < QCAD::MAX_FIELDNAMES_IN_INTEGRAL; i++) {
    fieldName = plist->get<std::string>(Albany::strint("Field Name",i),"");
    if(fieldName.length() > 0) {
      fieldNames.push_back(fieldName);
      conjugateFieldFlag.push_back(plist->get<bool>(Albany::strint("Conjugate Field",i),false));

      fieldName = plist->get<std::string>(Albany::strint("Field Name Im",i),"");
      fieldNames_Imag.push_back(fieldName);
      if(fieldName.length() > 0) fieldIsComplex.push_back(true);
      else fieldIsComplex.push_back(false);
    }
    else break;
  }
  bReturnImagPart = plist->get<bool>("Return Imaginary Part",false);
  
  std::string integrandLinLengthUnit; // linear length unit of integrand (e.g. "cm" for integrand in cm^-3)
  integrandLinLengthUnit = plist->get<std::string>("Integrand Length Unit","cm");
  bPositiveOnly = plist->get<bool>("Positive Return Only",false);

  //! compute scaling factor based on number of dimensions and units
  double integrand_length_unit_in_m;
  if( integrandLinLengthUnit == "m" )       integrand_length_unit_in_m = 1.0;
  else if( integrandLinLengthUnit == "cm" ) integrand_length_unit_in_m = 1e-2;
  else if( integrandLinLengthUnit == "um" ) integrand_length_unit_in_m = 1e-6;
  else if( integrandLinLengthUnit == "nm" ) integrand_length_unit_in_m = 1e-9;
  else if( integrandLinLengthUnit == "mesh" ) integrand_length_unit_in_m = length_unit_in_m;
  else integrand_length_unit_in_m = length_unit_in_m;  // assume same unit as mesh (e.g. if unit string is blank)
  
  double X0 = length_unit_in_m / integrand_length_unit_in_m; // length scaling to get to integrand's lenght unit

  if (numDims == 1)       scaling = X0; 
  else if (numDims == 2)  scaling = X0*X0; 
  else if (numDims == 3)  scaling = X0*X0*X0; 
  else 
    TEUCHOS_TEST_FOR_EXCEPTION (true, Teuchos::Exceptions::InvalidParameter, std::endl 
				<< "Error! Invalid number of dimensions: " << numDims << std::endl);


  //! add dependent fields (all fields assumed scalar qp)
  std::vector<std::string>::const_iterator it;
  //for(it = fieldNames.begin(); it != fieldNames.end(); ++it) {
  for(std::size_t i=0; i<fieldNames.size(); i++) {
    PHX::MDField<ScalarT,Cell,QuadPoint> f(fieldNames[i], scalar_dl);
    fields.push_back(f); this->addDependentField(f.fieldTag());
//.........这里部分代码省略.........