C++ Mesh::BeginFace方法代码示例

//.........这里部分代码省略.........
	{ // Prepare geometrical data on the mesh.
		Mesh::GeomParam table;
		table[CENTROID]    = CELL | FACE; //Compute averaged center of mass
		table[NORMAL]      = FACE;        //Compute normals
		table[ORIENTATION] = FACE;        //Check and fix normal orientation
		table[MEASURE]     = CELL | FACE; //Compute volumes and areas
		//table[BARYCENTER]  = CELL | FACE; //Compute volumetric center of mass
		mesh->RemoveGeometricData(table);
		mesh->PrepareGeometricData(table); //Ask to precompute the data
	}

	printf("nodes: %d edges: %d faces: %d cells: %d\n", mesh->NumberOfNodes(), mesh->NumberOfEdges(), mesh->NumberOfFaces(), mesh->NumberOfCells());


	{
		Storage::bulk_array name = mesh->self()->BulkArray(mesh->CreateTag("PROBLEMNAME",DATA_BULK,MESH,NONE));
		name.replace(name.begin(),name.end(),problem_name.begin(),problem_name.end());
	}

	Tag bndcond = mesh->CreateTag("BOUNDARY_CONDITION",DATA_REAL,FACE,FACE,3);
	Tag velocity = mesh->CreateTag("VELOCITY",DATA_REAL,FACE,NONE,1);
	Tag reaction = mesh->CreateTag("REACTION",DATA_REAL,CELL,NONE,1);
	Tag tensor = mesh->CreateTag("PERM", DATA_REAL, CELL, NONE, 1);
	Tag force = mesh->CreateTag("FORCE",DATA_REAL,CELL,NONE,1);
	Tag vel3 = mesh->CreateTag("VELVEC",DATA_REAL,CELL,NONE,3);
	Tag refsol = mesh->CreateTag("REFERENCE_SOLUTION",DATA_REAL,CELL,NONE,1);
	Tag refflux = mesh->CreateTag("REFERENCE_FLUX",DATA_REAL,FACE,NONE,1);
	Tag zone = mesh->CreateTag("ZONE",DATA_INTEGER,CELL,NONE,1);
	int numinner = 0, numouter = 0;
	int numinnern = 0, numoutern = 0;
	const double eps = 1.0e-6;
	const double mu = 1.0e-3;
	const double velmult = 1;
	for(Mesh::iteratorFace it = mesh->BeginFace(); it != mesh->EndFace(); ++it) 
	{
		Storage::real cnt[3], nrm[3];
		it->Centroid(cnt);
		it->UnitNormal(nrm);
		Storage::real x = cnt[0];
		Storage::real y = cnt[1];
		Storage::real z = cnt[2];
		Storage::real r = sqrt(x*x+y*y);
		Storage::real theta = atan2(y,x);//+pi;
		Storage::real sol, diff, flux, velnrm, dsolx, dsoly;
		Storage::real vel[3] = {-y/(r*r)*velmult,x/(r*r)*velmult,0.0};

		if( theta < 0 ) theta = 2*pi + theta;
		
		velnrm = nrm[0]*vel[0] + nrm[1]*vel[1] + nrm[2]*vel[2];

		if( theta < pi )
		{
			diff = pi;
			sol = (theta-pi)*(theta-pi);
			dsolx = -2*y*(theta-pi)/(r*r);
			dsoly = 2*x*(theta-pi)/(r*r);
			flux = velnrm*sol - diff*(dsolx*nrm[0] + dsoly*nrm[1]);
		}
		else
		{
			diff = 0;
			sol = 3*pi*(theta-pi);
			flux = velnrm*sol;
		}
		
		it->Real(refflux) = flux;

int main(int argc,char ** argv)
{
	Solver::Initialize(&argc,&argv,""); // Initialize the solver and MPI activity
#if defined(USE_PARTITIONER)
	Partitioner::Initialize(&argc,&argv); // Initialize the partitioner activity
#endif
	if( argc > 1 )
	{
		TagReal phi;
		TagReal tag_F;
		TagRealArray tag_K;
		TagRealArray tag_BC;
		TagReal phi_ref;
		Mesh * m = new Mesh(); // Create an empty mesh
		double ttt = Timer();
		bool repartition = false;
		m->SetCommunicator(INMOST_MPI_COMM_WORLD); // Set the MPI communicator for the mesh
		if( m->GetProcessorRank() == 0 ) // If the current process is the master one
			std::cout << argv[0] << std::endl;

		if( m->isParallelFileFormat(argv[1]) )
		{
			m->Load(argv[1]); // Load mesh from the parallel file format
			repartition = true;
		}
		else
		{
			if( m->GetProcessorRank() == 0 )
				m->Load(argv[1]); // Load mesh from the serial file format
		}
		BARRIER;
		if( m->GetProcessorRank() == 0 ) std::cout << "Processors: " << m->GetProcessorsNumber() << std::endl;
		if( m->GetProcessorRank() == 0 ) std::cout << "Load(MPI_File): " << Timer()-ttt << std::endl;

		//~ double ttt2 = Timer();
		//~ Mesh t;
		//~ t.SetCommunicator(INMOST_MPI_COMM_WORLD);
		//~ t.SetParallelFileStrategy(0);
		//~ t.Load(argv[1]);
		//~ BARRIER
		//~ if( m->GetProcessorRank() == 0 ) std::cout << "Load(MPI_Scatter): " << Timer()-ttt2 << std::endl;

#if defined(USE_PARTITIONER)
		if (m->GetProcessorsNumber() > 1)
		{ // currently only non-distributed meshes are supported by Inner_RCM partitioner
			ttt = Timer();
			Partitioner * p = new Partitioner(m);
			p->SetMethod(Partitioner::INNER_KMEANS,Partitioner::Partition); // Specify the partitioner
			p->Evaluate(); // Compute the partitioner and store new processor ID in the mesh
			delete p;
			BARRIER;

			if( m->GetProcessorRank() == 0 ) std::cout << "Evaluate: " << Timer()-ttt << std::endl;

			ttt = Timer();
			m->Redistribute(); // Redistribute the mesh data
			m->ReorderEmpty(CELL|FACE|EDGE|NODE); // Clean the data after reordring
			BARRIER;

			if( m->GetProcessorRank() == 0 ) std::cout << "Redistribute: " << Timer()-ttt << std::endl;
		}
#endif

		ttt = Timer();
		phi = m->CreateTag("Solution",DATA_REAL,CELL,NONE,1); // Create a new tag for the solution phi
		
		bool makerefsol = true;
		
		if( m->HaveTag("PERM" ) )
		{
			tag_K = m->GetTag("PERM");
			makerefsol = false;
			std::cout << "Permeability from grid" << std::endl;
		}
		else
		{
			std::cout << "Set perm" << std::endl;
			tag_K = m->CreateTag("PERM",DATA_REAL,CELL,NONE,1); // Create a new tag for K tensor
			for( Mesh::iteratorCell cell = m->BeginCell(); cell != m->EndCell(); ++cell ) // Loop over mesh cells
				tag_K[*cell][0] = 1.0; // Store the tensor K value into the tag
		}
		
		
		
		if( m->HaveTag("BOUNDARY_CONDITION") )
		{
			tag_BC = m->GetTag("BOUNDARY_CONDITION");
			makerefsol = false;
			std::cout << "Boundary conditions from grid" << std::endl;
		}
		else
		{
			std::cout << "Set boundary conditions" << std::endl;
			double x[3];
			tag_BC = m->CreateTag("BOUNDARY_CONDITION",DATA_REAL,FACE,FACE,3);
			for( Mesh::iteratorFace face = m->BeginFace(); face != m->EndFace(); ++face )
				if( face->Boundary() && !(face->GetStatus() == Element::Ghost) )
				{
					face->Centroid(x);
					tag_BC[*face][0] = 1; //dirichlet
//.........这里部分代码省略.........

int main(int argc, char ** argv)
{
	double nx = 2.0/7.0, ny = 6.0/7.0, nz = 3.0/7.0;
	double px = 0.5, py = 0.5, pz = 0.5;

	if( argc < 2 )
	{
		std::cout << "Usage: " << argv[0] << " mesh [mesh_out=grid.pmf] [nx=0] [ny=0] [nz=1] [px=0.5] [py=0.5] [pz=0.5]" << std::endl;
		return -1;
	}

	std::string grid_out = "grid.pmf";

	if( argc > 2 ) grid_out = std::string(argv[2]);
	if( argc > 3 ) nx = atof(argv[3]);
	if( argc > 4 ) ny = atof(argv[4]);
	if( argc > 5 ) nz = atof(argv[5]);
	if( argc > 6 ) px = atof(argv[6]);
	if( argc > 7 ) py = atof(argv[7]);
	if( argc > 8 ) pz = atof(argv[8]);

	double d = nx*px+ny*py+nz*pz;

	Mesh m;
	m.Load(argv[1]);
	m.SetTopologyCheck(NEED_TEST_CLOSURE|PROHIBIT_MULTILINE|PROHIBIT_MULTIPOLYGON|GRID_CONFORMITY|DEGENERATE_EDGE|DEGENERATE_FACE|DEGENERATE_CELL | FACE_EDGES_ORDER);
	//m.RemTopologyCheck(THROW_EXCEPTION);
	Tag sliced = m.CreateTag("SLICED",DATA_BULK,FACE|EDGE|NODE,FACE|EDGE|NODE,1);

	std::cout << "Cells: " << m.NumberOfCells() << std::endl;
	std::cout << "Faces: " << m.NumberOfFaces() << std::endl;

	MarkerType slice = m.CreateMarker();
	int nslice = 0, nmark = 0;
	for(Mesh::iteratorEdge it = m.BeginEdge(); it != m.EndEdge(); ++it)
	{
		double p[3];
		Storage::real_array c0 = it->getBeg()->Coords();
		Storage::real_array c1 = it->getEnd()->Coords();
		double r0 = c0[0]*nx+c0[1]*ny+c0[2]*nz - d;
		double r1 = c1[0]*nx+c1[1]*ny+c1[2]*nz - d;
		//std::cout << "r0 " << r0 << " r1 " << r1 << std::endl;
		if( r0*r1 < -1.0e-12 )
		{
			p[0] = (r0*c1[0] - r1*c0[0])/(r0-r1);
			p[1] = (r0*c1[1] - r1*c0[1])/(r0-r1);
			p[2] = (r0*c1[2] - r1*c0[2])/(r0-r1);
			//std::cout << "p " << p[0] << " " << p[1] << " " << p[2] << std::endl;
			Node n = m.CreateNode(p);
			n.Bulk(sliced) = 1;
			n.SetMarker(slice);
			bool was_sliced = it->HaveData(sliced) ? true : false;
			ElementArray<Edge> ret = Edge::SplitEdge(it->self(),ElementArray<Node>(&m,1,n.GetHandle()),0);
			if( was_sliced ) for(int q = 0; q < ret.size(); ++q) ret[q]->Bulk(sliced) = 1;
			nslice++;
		}
		else
		{
			if( fabs(r0) < 1.0e-6 )
			{
				it->getBeg()->SetMarker(slice);
				nmark++;
			}
			if( fabs(r1) < 1.0e-6 )
			{
				it->getEnd()->SetMarker(slice);
				nmark++;
			}
		}
	}

	std::cout << "sliced edges: " << nslice << " marked nodes: " << nmark << std::endl;

	if( !Element::CheckConnectivity(&m) )
		std::cout << "Connectivity is broken" << std::endl;
	
	nslice = 0;
	for(Mesh::iteratorFace it = m.BeginFace(); it != m.EndFace(); ++it)
	{
		ElementArray<Node> nodes = it->getNodes(slice); //those nodes should be ordered so that each pair forms an edge
		if( nodes.size() > 1 ) // if there is 1, then only one vertex touches the plane
		{
			//if there is more then two, then original face is non-convex
			if( nodes.size() > 2 ) std::cout << "Looks like face " << it->LocalID() << " is nonconvex" << std::endl;
			else
			{
				Edge e = m.CreateEdge(nodes).first;
				e.Bulk(sliced) = 1;
				e.SetMarker(slice);
				bool was_sliced = it->HaveData(sliced) ? true : false;
				ElementArray<Face> ret = Face::SplitFace(it->self(),ElementArray<Edge>(&m,1,e.GetHandle()),0);
				if( was_sliced ) for(int q = 0; q < ret.size(); ++q) ret[q]->Bulk(sliced) = 1;
				nslice++;
			}
		}
		//else std::cout << "Only one adjacent slice node, face " << it->LocalID() << std::endl;
	}

	nmark = 0;
	for(Mesh::iteratorEdge it = m.BeginEdge(); it != m.EndEdge(); ++it)
//.........这里部分代码省略.........