C++ GetPot::search方法代码示例

  std::string Runner::check_and_get_inputfile(int argc, char* argv[], GetPot & command_line)
  {
    if( argc < 2 )
      {
        std::stringstream error_msg;
        error_msg << "ERROR: Found only 1 command line argument, but was expecting an inputfile name!"
                  << std::endl
                  << "       Please specify the name of the input file on the command line as the first" << std::endl
                  << "       command line argument or using the '--input <filename>' option." << std::endl;
        libmesh_error_msg(error_msg.str());
      }

    std::string inputfile_name;
    if( command_line.search("--input") )
      inputfile_name = command_line.next(std::string("DIE!"));
    else
      inputfile_name = argv[1];

    std::ifstream i(inputfile_name.c_str());
    if (!i)
      {
        std::string error_msg = "Error: Could not read from input file "+inputfile_name+"!\n";
        libmesh_error_msg(error_msg);
      }

    return inputfile_name;
  }

T assert_argument (GetPot &cl,
                   const std::string &argname,
                   const char        *progname,
                   const T&          defaultarg)
{
  if(!cl.search(argname))
    {
      libMesh::err << ("No " + argname + " argument found!") << std::endl;
      usage_error(progname);
    }
  return cl.next(defaultarg);
}

int main (int argc, char *argv[])
	{
	const unsigned int ImageDimension = 2;
	    GetPot cl (argc, const_cast<char**>(argv));
	    if( cl.size() == 1 || cl.search (2,"--help","-h") )
	    {
	        std::cout << "Not Enough Arguments" << std::endl;
	        std::cout << "Generate the Gradient Table" << std::endl;
	        std::cout << "Usage:  return -1" << std::endl;
	    }

	    const string image_n = cl.follow("NoFile",1, "-i");
	    const string out_n   = cl.follow("NoFile",1, "-o");

	    typedef itk::DiffusionTensor3D<float> DiffusionTensorType;
	    typedef itk::Image<DiffusionTensorType, 3> ImageType;
	    typedef itk::ImageFileReader<ImageType> ReaderType;
	    ReaderType::Pointer reader = ReaderType::New();

	    reader->SetFileName(image_n);
	    reader->Update();

	    ImageType::Pointer image = reader->GetOutput();

	    typedef itk::ImageRegionIterator<ImageType> TensorIterator;
	    TensorIterator itImg(image, image->GetLargestPossibleRegion());

	    std::ofstream  file;
	    file.open(out_n);

	    for(itImg.GoToBegin(); !itImg.IsAtEnd(); ++itImg)
	    {
	    	file << itImg.Get() << std::endl;
	    }

	    file.close();



		return 0;
	}

int main (int argc, char **argv)
{

  GetPot cl (argc, argv);

  if (cl.search (2, "-h", "--help"))
    {
      std::cerr << help_text << std::endl;
      return 0;
    }
  const double a = cl.follow (0.0, "-a");
  const double b = cl.follow (1.0, "-b");
  const unsigned int nnodes = cl.follow (100, 2, "-n", "--nnodes");
  const std::string diffusion = cl.follow ("1.0", 2, "-d", "--diffusion");
  const std::string forcing = cl.follow ("1.0", 2, "-f", "--forcing");
  
  std::unique_ptr<mesh<double>> m (new mesh<double> (a, b, nnodes));
  fem_1d<double> problem (std::move (m));

  coeff<double> a_coeff (diffusion);
  problem.set_diffusion_coefficient (a_coeff);

  coeff<double> f_coeff (forcing);
  problem.set_source_coefficient (f_coeff);

  problem.assemble ();
  
  problem.set_dirichlet (fem_1d<double>::left_boundary, 0.0);
  problem.set_dirichlet (fem_1d<double>::right_boundary, 0.0);

  problem.solve ();
  
  for (unsigned int ii = 0; ii < nnodes; ++ii)
    std::cout << problem.m->nodes[ii] << " "
              << problem.result ()(ii, 0)
              << std:: endl;
      
  return 0;
};

void get_input_params(GetPot& field_parser,
		libmesh_assemble_input_params& input_params) {

	// Set mesh files
	if (field_parser.search(1, "Mesh")) {
		input_params.mesh_file = field_parser.next(
				input_params.mesh_file);
	} else {
		homemade_error_msg("Missing the system mesh file!");
	}

	// Set constant parameters
	if ( field_parser.search(1, "PhysicalParameters") )
	{
		input_params.physical_params_file = field_parser.next(input_params.physical_params_file);
	}
	else
	{
		homemade_error_msg("Missing the physical parameters file!");
	}

	// Set weight function
	std::string sys_type;
	if ( field_parser.search(1, "SystemType") )
	{
		sys_type = field_parser.next(sys_type);
		if(sys_type == "Macro" || sys_type == "MACRO" || sys_type == "macro")
			input_params.system_type = WeightFunctionSystemType::MACRO;
		else if(sys_type == "Micro" || sys_type == "MICRO" || sys_type == "micro")
			input_params.system_type = WeightFunctionSystemType::MICRO;
		else if(sys_type == "NoWeight" || sys_type == "NOWEIGHT" || sys_type == "noweight")
		{
			input_params.system_type = WeightFunctionSystemType::NO_WEIGHT;
			std::cout << " >> Warning: Will not use the weight parameters!" << std::endl;
		}
		else
			homemade_error_msg("Invalid system type (must be either Macro, Micro or NoWeight)!");
	}
	else
	{
		homemade_error_msg("Missing the system type (must be either Macro, Micro or NoWeight)!");
	}

	if ( field_parser.search(1, "MeshWeight") )
	{
		input_params.mesh_weight_file = field_parser.next(input_params.mesh_weight_file);
	}
	else
	{
		homemade_error_msg("Missing the weight mesh file!");
	}

	if( field_parser.search(1, "WeightIndexes") )
	{
		input_params.weight_domain_idx_file = field_parser.next(input_params.weight_domain_idx_file);
	}
	else
	{
		homemade_error_msg("Missing the weight value file!");
	}

	// Output
	if (field_parser.search(2, "--output", "OutputBase"))
	{
		input_params.output_base = field_parser.next(
			input_params.output_base);
	} else {
		input_params.output_base = "test_system";
	}

	if (field_parser.search(1, "ExportRBVectors")) {
		input_params.bCalculateRBVectors = true;
	} else {
		input_params.bCalculateRBVectors = false;
	}
};

int main (int argc, char **argv)
{

    GetPot cl (argc, argv);

    if (cl.search (2, "-h", "--help"))
    {
        std::cerr << help_text << std::endl;
        return 0;
    }
    const double a = cl.follow (0.0, "-a");
    const double b = cl.follow (1.0, "-b");
    const unsigned int nnodes = cl.follow (100, 2, "-n", "-nnodes");

    mesh m (a, b, nnodes);
    Eigen::MatrixXd A(nnodes, nnodes);

    Eigen::Matrix2d mloc;
    mloc << 0, 0, 0, 0;
    for (unsigned int iel = 0; iel < m.nels; ++iel)
    {

        mloc << 0, 0, 0, 0;
        for (unsigned int inode = 0; inode < 2; ++inode)
        {
            double igrad = (inode == 0 ? 1.0 / m.h : -1.0 / m.h);
            for (unsigned int jnode = 0; jnode < 2; ++jnode)
            {
                double jgrad =  (jnode == 0 ? 1.0 / m.h : -1.0 / m.h);
                mloc(inode,jnode) = igrad * jgrad * m.h;
                A(m.elements[iel][inode],m.elements[iel][jnode]) +=
                    mloc(inode,jnode);
            }
        }
    }

    Eigen::VectorXd f(nnodes);
    Eigen::Vector2d vloc;

    for (unsigned int iel = 0; iel < m.nels; ++iel)
    {
        vloc << 0, 0;

        for (unsigned int inode = 0; inode < 2; ++inode)
        {
            vloc(inode) = m.h / 2.0;
            f(m.elements[iel][inode]) += vloc(inode);
        }
    }

    f(0) = 0;
    f(nnodes - 1) = 0;

    A(0,0) = 1.0;
    A(nnodes-1,nnodes-1) = 1.0;
    for (unsigned int ii = 1; ii < nnodes; ++ii)
    {
        A(0, ii) = 0.0;
        A(nnodes-1, nnodes-1-ii) = 0.0;
    }


    Eigen::VectorXd uh = A.partialPivLu ().solve (f);

    for (unsigned int ii = 0; ii < nnodes; ++ii)
        std::cout << uh(ii, 0) << std:: endl;

    return 0;
};

int main (int argc, char **argv)
{

  GetPot cl (argc, argv);

  if (cl.search (2, "-h", "--help"))
    {
      std::cerr << help_text << std::endl;
      return 0;
    }
  const double a = cl.follow (0.0, "-a");
  const double b = cl.follow (1.0, "-b");
  const unsigned int nnodes = cl.follow (100, 2, "-n", "-nnodes");

  mesh m (a, b, nnodes);
  matrix A(nnodes);

  matrix mloc(2);
  for (unsigned int iel = 0; iel < m.nels; ++iel)
    {

      std::fill (mloc.get_data (),
                 mloc.get_data () + 4,
                 0.0);
      
      for (unsigned int inode = 0; inode < 2; ++inode)
        {
          double igrad = (inode == 0 ? 1.0 / m.h : -1.0 / m.h);
          for (unsigned int jnode = 0; jnode < 2; ++jnode)
            {
              double jgrad =  (jnode == 0 ? 1.0 / m.h : -1.0 / m.h);
              mloc(inode,jnode) = igrad * jgrad * m.h;
              A(m.elements[iel][inode],m.elements[iel][jnode]) += 
                mloc(inode,jnode);
            }
        }
    }

  matrix f(nnodes, 1);
  matrix vloc(2, 1);
  
  for (unsigned int iel = 0; iel < m.nels; ++iel)
    {
      std::fill (mloc.get_data (),
                 mloc.get_data () + 2,
                 0.0);
      
      for (unsigned int inode = 0; inode < 2; ++inode)
        {
          vloc(inode, 0) = m.h / 2.0;
          f(m.elements[iel][inode], 0) += vloc(inode, 0);
        }
    }

  f(0, 0) = 0;
  f(nnodes - 1, 0) = 0;

  A(0,0) = 1.0;
  A(nnodes-1,nnodes-1) = 1.0;
  for (unsigned int ii = 1; ii < nnodes; ++ii)
    {
      A(0, ii) = 0.0;
      A(nnodes-1, nnodes-1-ii) = 0.0;
    }

  
  matrix uh(f);
  A.solve (uh);

  for (unsigned int ii = 0; ii < nnodes; ++ii)
    std::cout << uh(ii, 0) << std:: endl;
      
  return 0;
};

void Server::parseConfig(GetPot &config)
{
  // this variable is needed for setting up the partial loading part when the features
  // are loaded for the very first time at fire start up
  string partialLoadingString="empty";

  if(config.search(2,"-F","--filter"))
  {
    string unparsed = config.follow("empty",2,"-F","--filter");
    if (parseFilter(unparsed.c_str()))
    {
      retriever_.setFilterApply(true);
      DBG(10) << "Using filter " << unparsed << " for retrieval" << endl;
    }
    else
    {
      retriever_.setFilterApply(false);
      retriever_.clearFilter();
      DBG(10) << "Error in filter String " << unparsed << endl;
      DBG(10) << "Aborting filtering and using normal retrieval settings" << endl;
    }
  }

  if(config.search(2,"-q","--queryCombiner")) {
    retriever_.setQueryCombiner(config.follow("adding:POS=1.0:NEG=0.833",2,"-q","--queryCombiner"));
  }

  if(config.search("--reRanker")) {
    retriever_.setReranking(config.follow("cluster:CONS=100:RR=20:CLUSTERS=5","--reRanker"));
  } 
  if(config.search(2,"-U","--defdontload"))
  {
    partialLoadingString="default";
    retriever_.setPartialLoadingApply(true);
  }

  if(config.search(2,"-u","--dontload"))
  {
    string unparsed=config.follow("empty",2,"-u","--dontload");
    if(unparsed!="empty")
    {
      // the parsing is not possible here, because the amount of different features
      // is not known yet
      // but a sanity check is possible
      const char* str = unparsed.c_str();
      bool errformat = false;
      // note it is not checked if there are two or more : present in a row
      while(*str != '\0' && !errformat)
      {
        if(!isdigit(*str) && *str!=':')
        {
          errformat=true;
          DBG(10) << "Error in format string " << unparsed << " following -u/--dontload option" << endl;
          DBG(10) << "format string must start with a digit; aborting partial feature loading" << endl;
        }
        str++;
      }
      if(!errformat && retriever_.setPartialLoadingApply(true))
      {
        partialLoadingString=unparsed;
      }
      else
      {
        ERR << "-u/--dontload was used without -F/--filter option or there are illegal characters in the format string" << endl;
        ERR << "ignoring -u/--dontload option, loading all features" << endl;
      }
    }
    else
    {
      retriever_.setPartialLoadingApply(false);
      ERR << "Error after -u/--dontload option: No featureindexsequence given" << endl;
      //DBG(10) << "Aborting partial loading for filtered retrieval and loading all features" << endl;
    }
  }

  if(config.search(2,"-f","--filelist"))
  {
    string filelistname=config.follow("list.txt",2,"-f","--filelist");
    string result=retriever_.filelist(filelistname,partialLoadingString);
    if(result=="filelist FAILURE")
    {
      ERR << "Error in filelist. Please correct it and start again." << endl;
      exit(1);
    }
    DBG(10) << result << endl;
  }
  if(config.search(2,"-R","--relevancefile"))
  {
    relevanceFile_=config.follow("",2,"-R","--relevancefile");
  }

  retriever_.setScoring("linear");

  
  log_=LogFile(config.follow("",2,"-l","--logfile"));

  if(config.search(2,"-D","--defaultdists"))
  {
    for(uint i=0;i<retriever_.numberOfSuffices();++i)
    {
//.........这里部分代码省略.........

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

    GetPot cl (argc, const_cast<char**>(argv));
    if( cl.size() == 1 || cl.search (2,"--help","-h") )
    {
        std::cout << "Not Enough Arguments" << std::endl;
        std::cout << "Scales the tensors with a scalar factor" << std::endl;
        std::cout << "Usage: -trueB0 <true B0> -m <MaskImage> -true <True Tensors> -f <flag for extended gradient> -t <initial tensor estimate> -g <gradient> -o <Output File> -s <Sigma> -nm <Noise Model> -Sim <intelligent COnvergence>" << std::endl;
      return -1;
    }

    const string image_n =cl.follow("NoFile", 1, "-i");
    const string out_n = cl.follow("NoFile", 1, "-o");

    const int idx = cl.follow(0, 1, "-ix");
    const int idy = cl.follow(0, 1, "-iy");
    const int idz = cl.follow(0, 1, "-iz");

    const int sx = cl.follow(0, 1, "-sx");
    const int sy = cl.follow(0, 1, "-sy");
    const int sz = cl.follow(0, 1, "-sz");



    typedef itk::Image< float, 3> ImageType;
    typedef itk::ImageFileReader<ImageType> ReaderType;
    ReaderType::Pointer reader = ReaderType::New();

    reader->SetFileName(image_n);
    reader->Update();

    ImageType::Pointer image = reader->GetOutput();

    ImageType::SizeType size;
    ImageType::IndexType id;

    id[0] = idx; id[1] = idy; id[2] = idz;

    size[0] = sx; size[1] = sy; size[2] =sz;

    ImageType::RegionType region(id, size);

    typedef itk::ExtractImageFilter<ImageType, ImageType> FilterType;
    FilterType::Pointer filter = FilterType::New();

    filter->SetExtractionRegion(region);
    filter->SetInput(image);
    filter->SetDirectionCollapseToIdentity();

    filter->Update();

    ImageType::Pointer exImage = filter ->GetOutput();

    typedef itk::ImageFileWriter<ImageType> WriterType;
    WriterType::Pointer writer = WriterType::New();

    writer->SetFileName(out_n);
    writer->SetInput(exImage);
    writer->Update();






    return 0;
}

int main(int argc, char *argv[])
{
    GetPot cl (argc, const_cast<char**>(argv));
    if( cl.size() == 1 || cl.search (2,"--help","-h") )
	    {
	        std::cout << "Not Enough Arguments" << std::endl;
	        std::cout << "Generate the Gradient Table" << std::endl;
	        std::cout << "Usage:  return -1" << std::endl;
	    }

    // One idea is to apply the rotation matrix
    const string Linear_trans_n = cl.follow("NoFile",1,"-t");
    const string ref_Image_n = cl.follow("NoFile",1, "-r");
    const string out_n = cl.follow("NoFile",1,"-o");
    const string B0_n =cl.follow("NoFile",1,"-B0");
    const string T1_n =cl.follow("NoFile",1,"-T1");

    typedef itk::Image<float, 3> ImageType;
    typedef itk::ImageFileReader<ImageType> ImageReaderType;
    ImageReaderType::Pointer imageReaderB0 = ImageReaderType::New();
    imageReaderB0->SetFileName(B0_n);
    imageReaderB0->Update();

    std::cout << imageReaderB0->GetOutput()->GetDirection() << std::endl;


    ImageReaderType::Pointer imageReaderT1 = ImageReaderType::New();

    imageReaderT1->SetFileName(B0_n);
        imageReaderT1->Update();


        std::cout << imageReaderT1->GetOutput()->GetDirection() << std::endl;


//    typedef itk::TransformFileReader TransformFileReaderType;
//    typedef TransformFileReaderType::TransformListType TransformListType;
//    typedef itk::TransformBase TransformBaseType;
//    typedef itk::AffineTransform<double, 3> AffineTransformType;
//
//    typedef itk::Image<float, 3> ImageType;
//    typedef itk::ImageFileReader<ImageType> ImageReaderType;
//    ImageReaderType::Pointer imageReader = ImageReaderType::New();
//
//    imageReader->SetFileName(ref_Image_n);
//    imageReader->Update();
//    ImageType::Pointer refImage = imageReader->GetOutput();
//
//
//
//	TransformFileReaderType::Pointer readerTransform = TransformFileReaderType::New();
//	readerTransform->SetFileName(Linear_trans_n);
//	readerTransform -> Update();
//	TransformListType *list = readerTransform->GetTransformList();
//	TransformBaseType * transform = list->front().GetPointer();
//	TransformBaseType::ParametersType parameters = transform->GetParameters();
//	AffineTransformType::Pointer transform_fwd = AffineTransformType::New();
//	transform_fwd->SetParameters(parameters);
//
//	std::cout << transform_fwd->GetParameters() << std::endl;

//    typedef itk::Vector< float, 3 >          VectorPixelType;
//    typedef itk::Image< VectorPixelType, 3 > DisplacementFieldImageType;
//
//
//    typedef itk::TransformToDisplacementFieldFilter<DisplacementFieldImageType, double> DisplacementFieldGeneratorType;
//    DisplacementFieldGeneratorType::Pointer dispfieldGenerator = DisplacementFieldGeneratorType::New();
//
//    dispfieldGenerator->UseReferenceImageOn();
//    dispfieldGenerator->SetReferenceImage( refImage );
//    dispfieldGenerator->SetTransform( transform_fwd );
//    dispfieldGenerator->Update();
//    DisplacementFieldImageType::Pointer dispField = dispfieldGenerator->GetOutput();
//
//    typedef itk::ImageFileWriter<DisplacementFieldImageType> DispFieldWriterType;
//    DispFieldWriterType::Pointer writer = DispFieldWriterType::New();
//    writer->SetFileName(out_n);
//    writer->SetInput(dispField);
//    writer->Update();

    return 0;
}

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

    GetPot cl (argc, const_cast<char**>(argv));
    if( cl.size() == 1 || cl.search (2,"--help","-h") )
    {
        std::cout << "Not Enough Arguments" << std::endl;
        std::cout << "Scales the tensors with a scalar factor" << std::endl;
        std::cout << "Usage: -trueB0 <true B0> -m <MaskImage> -true <True Tensors> -f <flag for extended gradient> -t <initial tensor estimate> -g <gradient> -o <Output File> -s <Sigma> -nm <Noise Model> -Sim <intelligent COnvergence>" << std::endl;
      return -1;
    }




     const string file_g_n = cl.follow("NoFile",1, "-g");
     const string fileIn  = cl.follow("NoFile",1,"-i");
     const string fileIn_HR  = cl.follow("NoFile",1,"-iHR");
     const string B0_n     =  cl.follow("NoFile", 1, "-B0");
     const string mask_LR_n   = cl.follow("NoFile",1, "-mLR");
     const int numOfIter   = cl.follow(1,1, "-n");
     const float kappa_L   = cl.follow(0.05, 1, "-k");
     const float lambda_L  = cl.follow(0.25, 1, "-lamb_L");
     const string B0Image_HR_n = cl.follow("NoFile", 1, "-B0HR");
     const string dispField_n = cl.follow("NoFile",1,"-d");
     const string T1Image_n   = cl.follow("NoFile",1,"-T1");
     const string mask_HR_n = cl.follow("NoFile", 1, "-mHR");
    // Usual Typedefs
    typedef float RealType;
    const int ImageDim  =3;

    typedef itk::Image<RealType, ImageDim> ScalarImageType;
    typedef itk::Vector<double, ImageDim> VectorType;
    typedef itk::Image<VectorType, ImageDim> DeformationFieldType;
    typedef itk::Image<VectorType, ImageDim> VectorImageType;

    typedef itk::ImageFileReader<ScalarImageType> ScalarImageReaderType;
    typedef itk::ImageFileWriter<ScalarImageType> ScalarImageWriterType;

    //Read T1 image
	ScalarImageReaderType::Pointer scalarReader = ScalarImageReaderType::New();
	scalarReader->SetFileName(T1Image_n.c_str());
	scalarReader->Update();
	
	ScalarImageType::Pointer T1_image = scalarReader->GetOutput();

    // Read LR ImageList
	typedef std::vector<ScalarImageType::Pointer> ImageListType;
	ImageListType DWIList;
	
	std::ifstream file(fileIn.c_str());
        int numOfImages = 0;
        file >> numOfImages;	
    
	for (int i=0; i < numOfImages  ; i++) // change of numOfImages
         {
             char filename[256];
             file >> filename;
             ScalarImageReaderType::Pointer myReader=ScalarImageReaderType::New();
             myReader->SetFileName(filename);
             std::cout << "Reading.." << filename << std::endl; // add a try catch block
             myReader->Update();
             DWIList.push_back( myReader->GetOutput() ); //using push back to create a stack of diffusion images
         }

	// Read deformation field
	typedef itk::ImageFileReader<DeformationFieldType> DeformationFieldReaderType;
	DeformationFieldReaderType::Pointer deformationFieldReader = DeformationFieldReaderType::New();

	deformationFieldReader->SetFileName(dispField_n.c_str());
	deformationFieldReader->Update();
	
	DeformationFieldType::Pointer defField = deformationFieldReader->GetOutput();
	

	// Read Mask Image Spatial
	typedef itk::ImageMaskSpatialObject<ImageDim> MaskSpatialObjectType;
	typedef MaskSpatialObjectType::ImageType MaskSpatialImageType;
	typedef itk::ImageFileReader<MaskSpatialImageType> MaskSpatialImageReaderType;	

	MaskSpatialImageReaderType::Pointer spatialReader = MaskSpatialImageReaderType::New();
	spatialReader->SetFileName(mask_LR_n.c_str());
	spatialReader->Update();
	
	MaskSpatialImageType::Pointer maskSpatialImage_LR = spatialReader->GetOutput();

	//Read Mask Image Normal
	ScalarImageReaderType::Pointer maskImageReader = ScalarImageReaderType::New();
	maskImageReader->SetFileName(mask_LR_n.c_str());
	maskImageReader->Update();
	
	ScalarImageType::Pointer maskImage_LR = maskImageReader->GetOutput();

		
	// Resample diffusion Images	
/*	typedef itk::WarpImageFilter<ScalarImageType, ScalarImageType, DeformationFieldType> WarpImageFilterType;
	typedef itk::ImageFileWriter<ScalarImageType> ScalarImageWriterType;
	
	

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

int main (int argc, char **argv)
{

  GetPot cl (argc, argv);

  if (cl.search (2, "-h", "--help"))
    {
      std::cerr << help_text << std::endl;
      return 0;
    }

  const double a =
    cl.follow (0.0, "-a");

  const double b =
    cl.follow (1.0, "-b");

  const unsigned int nnodes =
    cl.follow (100, 2, "-n", "--nnodes");

  const std::string diffusion =
    cl.follow ("1.0", 2, "-d", "--diffusion");

  const std::string forcing =
    cl.follow ("1.0", 2, "-f", "--forcing");
  
  coeff f_coeff (forcing);
  coeff a_coeff (diffusion);

  const std::string quadrature =
    cl.follow ("trapezoidal.so",
               2, "-q", "--quadrature-rule");

  std::function <void ()> r_r;
  
  void * handle = dlopen (quadrature.c_str (), RTLD_NOW);
  if (! handle)
    {
      std::cerr << "fem1d: cannot load dynamic object!"
                << std::endl;
      std::cerr << dlerror ()
                << std::endl;
      return (-1);
    }

  void * sym = dlsym (handle, "register_rules");
  if (! sym)
    {
      std::cerr << "fem1d: cannot load symbol!"
                << std::endl;
      std::cerr << dlerror ()
                << std::endl;
      return (-1);
    }


  r_r = reinterpret_cast <void (*) ()> (sym);
  r_r ();

  auto & the_factory = quadrature_factory::instance (); 
  auto integrate = the_factory.create ("trapezoidal");
  if (! integrate)
    {
      std::cerr << "rule name unknown" << std::endl;
      return (-1);
    }
    
  mesh m (a, b, nnodes);
  Eigen::SparseMatrix<double> A(nnodes, nnodes);
  
  Eigen::Matrix2d mloc;
  mloc << 0, 0, 0, 0;
  for (unsigned int iel = 0; iel < m.nels; ++iel)
    {

      mloc << 0, 0, 0, 0;      
      for (unsigned int inode = 0; inode < 2; ++inode)
        {

          auto igrad = [=, &m] (double x) -> double
            {
              return basisfungrad
              (x, m.nodes[m.elements[iel][0]],
               m.nodes[m.elements[iel][1]],
               inode == 0 ? 1.0 : 0.0, 
               inode == 1 ? 1.0 : 0.0);
            };

          
          for (unsigned int jnode = 0; jnode < 2; ++jnode)
            {

              auto jgrad = [=, &m] (double x) -> double
                {
                  return basisfungrad
                  (x, m.nodes[m.elements[iel][0]],
                   m.nodes[m.elements[iel][1]],
                   jnode == 0 ? 1.0 : 0.0, 
                   jnode == 1 ? 1.0 : 0.0);
                };
//.........这里部分代码省略.........

int main (int argc, char *argv[])
	{
	const unsigned int ImageDimension = 2;
	    GetPot cl (argc, const_cast<char**>(argv));
	    if( cl.size() == 1 || cl.search (2,"--help","-h") )
	    {
	        std::cout << "Not Enough Arguments" << std::endl;
	        std::cout << "Generate the Gradient Table" << std::endl;
	        std::cout << "Usage:  return -1" << std::endl;
	    }

	   const string image_n = cl.follow("NoFile",1, "-i");
	   const string mask_n = cl.follow("NoFile", 1, "-m");	 
	   const string out_n   = cl.follow("NoFile",1, "-o");

	    typedef itk::DiffusionTensor3D<float> DiffusionTensorType;
	    typedef itk::Image<DiffusionTensorType, 3> TensorImageType;
	    typedef itk::ImageFileReader<TensorImageType> TensorReaderType;
	    TensorReaderType::Pointer reader = TensorReaderType::New();

	    reader->SetFileName(image_n.c_str());
	    reader->Update();

	    TensorImageType::Pointer image = reader->GetOutput();


	    typedef itk::Image<float, 3> ScalarImageType;
	    typedef itk::ImageFileReader<ScalarImageType> ScalarReaderType;

	    ScalarReaderType::Pointer scalarReader = ScalarReaderType::New();
	    scalarReader->SetFileName(mask_n.c_str());
	    scalarReader->Update();
	    ScalarImageType::Pointer maskImage = scalarReader->GetOutput();
	
	    typedef itk::ImageRegionIterator<TensorImageType> TensorIterator;

	    typedef itk::ImageRegionIterator<ScalarImageType> ScalarIterator;
	    ScalarIterator itMask(maskImage, maskImage->GetLargestPossibleRegion());
	

	    TensorUtilities  utilsTensor;
	    TensorImageType::Pointer logTensorImage = utilsTensor.LogTensorImageFilter(image, maskImage);


	    typedef itk::ImageFileWriter<TensorImageType> TensorImageWriterType;
	    TensorImageWriterType::Pointer tensorImageWriter = TensorImageWriterType::New();

            tensorImageWriter->SetFileName("LogTensorImage_stupid.nii.gz");
	    tensorImageWriter->SetInput(logTensorImage);
	    tensorImageWriter->Update();
	
	    std::ofstream  file;
	    file.open(out_n);

	 	ScalarImageType::Pointer TraceImage = ScalarImageType::New();
	    CopyImage	cpImage;
	    cpImage.CopyScalarImage(maskImage, TraceImage);

	    ScalarIterator itTr(TraceImage, TraceImage->GetLargestPossibleRegion());



	    TensorIterator itImg(logTensorImage, logTensorImage->GetLargestPossibleRegion());
	    for(itImg.GoToBegin(), itMask.GoToBegin(), itTr.GoToBegin(); !itTr.IsAtEnd(), !itImg.IsAtEnd(), !itMask.IsAtEnd(); ++itTr, ++itImg, ++itMask)
	    {
	    		file << itImg.Get().GetTrace() << std::endl;
	    		itTr.Set(itImg.Get().GetTrace()) ;
		

	   }


		typedef itk::ImageFileWriter<ScalarImageType> WriterType;
		WriterType::Pointer writer = WriterType::New();
	
		writer->SetFileName("LogTense_Trace.nii.gz");
		writer->SetInput(TraceImage);
		writer->Update();


	    file.close();


		return 0;
	}

int main (int argc, char **argv)
{

  GetPot cl (argc, argv);

  if (cl.search (2, "-h", "--help"))
    {
      std::cerr << help_text << std::endl;
      return 0;
    }
  const double a = cl.follow (double (0.0), "-a");
  const double b = cl.follow (double (1.0), "-b");
  const unsigned int nnodes = cl.follow (100, 2, "-n", "--nnodes");
  const unsigned int nel = nnodes - 1;
  const std::string diffusion = cl.follow ("1.0", 2, "-d", "--diffusion");
  const std::string forcing = cl.follow ("1.0", 2, "-f", "--forcing");

  const double tol = 1e-6;
  const unsigned int maxit = 43;
  const unsigned int overlap = 3;
  const double L = b - a;
  const int mpi_size = 3;

  const double L_loc = L / double (mpi_size);
  const double h = L_loc / ceil (double(nel) / double(mpi_size));

  std::vector<double> a_loc (mpi_size);
  std::vector<double> b_loc (mpi_size);
  std::vector<unsigned int> nel_loc (mpi_size);
  std::vector<unsigned int> ndof_loc (mpi_size);
  std::vector<fem_1d<double>*> subproblems (mpi_size);
  
  coeff<double> a_coeff (diffusion);
  coeff<double> f_coeff (forcing);

  int mpi_rank;
  for (mpi_rank = 0; mpi_rank < mpi_size; ++mpi_rank)
    {
      a_loc[mpi_rank] = a + mpi_rank * L_loc;
      b_loc[mpi_rank] = a_loc[mpi_rank] + L_loc;
      nel_loc[mpi_rank] = ceil (double(nel) / double(mpi_size));
      if (mpi_rank > 0)
        {
          a_loc[mpi_rank] -= overlap * h;
          nel_loc[mpi_rank] += overlap;
        }
      if (mpi_rank < mpi_size - 1)
        {
          b_loc[mpi_rank] += overlap * h;
          nel_loc[mpi_rank] += overlap;
        }
      ndof_loc[mpi_rank] = nel_loc[mpi_rank] + 1;
      fem_1d<double>* tmp = new fem_1d<double>(new mesh<double>
                                               (a_loc[mpi_rank],
                                                b_loc[mpi_rank],
                                                ndof_loc[mpi_rank]));
      subproblems[mpi_rank] = tmp;
      
      subproblems[mpi_rank]->set_diffusion_coefficient (a_coeff);
      subproblems[mpi_rank]->set_source_coefficient (f_coeff);
      subproblems[mpi_rank]->assemble ();
      subproblems[mpi_rank]->set_dirichlet (fem_1d<double>::left_boundary, 0.0);
      subproblems[mpi_rank]->set_dirichlet (fem_1d<double>::right_boundary, 0.0);
      subproblems[mpi_rank]->solve ();
    };

  for (unsigned int it = 0; it < maxit; ++it)
    for (mpi_rank = 0; mpi_rank < mpi_size; ++mpi_rank)
      {
        if (mpi_rank > 0)
          subproblems[mpi_rank]->set_dirichlet
            (fem_1d<double>::left_boundary,
             subproblems[mpi_rank-1]->result ()
             [ndof_loc[mpi_rank-1]-1-2*overlap]);
        else
          subproblems[mpi_rank]->set_dirichlet
            (fem_1d<double>::left_boundary, 0.0);
        if (mpi_rank < mpi_size - 1)
          subproblems[mpi_rank]->set_dirichlet
            (fem_1d<double>::right_boundary,
             subproblems[mpi_rank+1]->result ()
             [2*overlap]);
        else
          subproblems[mpi_rank]->set_dirichlet
            (fem_1d<double>::right_boundary, 0.0);
        subproblems[mpi_rank]->solve ();
      }

  for (mpi_rank = 0; mpi_rank < mpi_size; ++mpi_rank)
    for (unsigned int ii = 0; ii < ndof_loc[mpi_rank]; ++ii)
      std::cout << subproblems[mpi_rank]->m->nodes[ii] << " "
                << subproblems[mpi_rank]->result ()(ii, 0)
                << std::endl;

  return 0;
};

Int
main ( Int argc, char** argv )
{
    //Setup main communicator
    boost::shared_ptr< Epetra_Comm > comm;

    //Setup MPI variables
    Int numberOfProcesses (1);
    Int rank (0);

#ifdef HAVE_MPI
    MPI_Init ( &argc, &argv );

    MPI_Comm_size ( MPI_COMM_WORLD, &numberOfProcesses );
    MPI_Comm_rank ( MPI_COMM_WORLD, &rank );
#endif

    if ( rank == 0 )
    {
        std::cout << std::endl;
        std::cout << "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$" << std::endl;
        std::cout << " THE ZERO DIMENSIONAL SOLVER IS AN ALPHA VERSION UNDER STRONG DEVELOPMENT" << std::endl;
        std::cout << "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$" << std::endl << std::endl;

        std::cout << "MPI Processes: " << numberOfProcesses << std::endl;
    }

#ifdef HAVE_MPI
    if ( numberOfProcesses > 1 )
    {
        if ( rank == 0 )
        {
            std::cout << "test_ZeroDimensional not enabled in parallel, failing gracefully." << std::endl;
            std::cout << "MPI Finalization" << std::endl;
        }
        MPI_Finalize();
        return EXIT_FAILURE;
    }
#endif

#ifdef EPETRA_MPI
    if ( rank == 0 )
    {
        std::cout << "MPI Epetra Initialization ... " << std::endl;
    }
    comm.reset ( new Epetra_MpiComm ( MPI_COMM_WORLD ) );
#else
    std::cout << "SERIAL Epetra Initialization ... " << std::endl;
    comm.reset ( new Epetra_SerialComm() );
#endif

    bool exitFlag = EXIT_SUCCESS;

#if ( defined(HAVE_NOX_THYRA) && defined(HAVE_TRILINOS_RYTHMOS) )
    // Command line parameters
    GetPot commandLine ( argc, argv );
    const bool check = commandLine.search ( 2, "-c", "--check" );
    string fileName  = commandLine.follow ( "data", 2, "-f", "--file" );

    // SetupData
    GetPot dataFile ( fileName  + ".dat" );

    std::string circuitDataFile = dataFile ( "0D_Model/CircuitDataFile", "./inputFile.dat" );
    BCInterface0D< ZeroDimensionalBCHandler, ZeroDimensionalData >  zeroDimensionalBC;
    zeroDimensionalBC.createHandler();
    zeroDimensionalBC.fillHandler ( circuitDataFile, "Files" );

    boost::shared_ptr< ZeroDimensionalData > zeroDimensionalData ( new ZeroDimensionalData );
    zeroDimensionalData->setup ( dataFile, zeroDimensionalBC.handler() );

    boost::shared_ptr< ZeroDimensionalSolver > zeroDimensionalSolver ( new ZeroDimensionalSolver ( zeroDimensionalData->unknownCounter(), comm, zeroDimensionalData->circuitData() ) );
    zeroDimensionalSolver->setup ( zeroDimensionalData->solverData() );

    zeroDimensionalData->showMe();

    // SetupModel
    zeroDimensionalData->dataTime()->setInitialTime (0);
    zeroDimensionalData->initializeSolution();

    // Create output folder
    if ( comm->MyPID() == 0 )
    {
        mkdir ( "output", 0777 );
    }

    // Save initial solution
    zeroDimensionalData->saveSolution();

    zeroDimensionalData->dataTime()->updateTime();
    zeroDimensionalData->dataTime()->setInitialTime (zeroDimensionalData->dataTime()->time() );

    // Definitions for the time loop
    LifeChrono chronoTotal;
    LifeChrono chronoSystem;
    LifeChrono chronoIteration;

    Int count = 0;
    chronoTotal.start();

    for ( ; zeroDimensionalData->dataTime()->canAdvance() ; zeroDimensionalData->dataTime()->updateTime(), ++count )
//.........这里部分代码省略.........