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

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


		const INMOST_DATA_INTEGER_TYPE nvf[24] = { 2, 3, 1, 0, 4, 5, 7, 6, 0, 1, 5, 4, 3, 2, 6, 7, 2, 0, 4, 6, 1, 3, 7, 5 };
		const INMOST_DATA_INTEGER_TYPE numnodes[6] = { 4, 4, 4, 4, 4, 4 };
		for (int i = 1; i < n; i++) 
		{
			for (int j = 1; j < n; j++) 
			{
				for (int k = 1; k < n; k++) 
				{

					ElementArray<Node> verts(mesh,8); 
					verts[0] = mesh->NodeByLocalID(V_ID(i - 1, j - 1, k - 1));
					verts[1] = mesh->NodeByLocalID(V_ID(i - 0, j - 1, k - 1));
					verts[2] = mesh->NodeByLocalID(V_ID(i - 1, j - 0, k - 1));
					verts[3] = mesh->NodeByLocalID(V_ID(i - 0, j - 0, k - 1));
					verts[4] = mesh->NodeByLocalID(V_ID(i - 1, j - 1, k - 0));
					verts[5] = mesh->NodeByLocalID(V_ID(i - 0, j - 1, k - 0));
					verts[6] = mesh->NodeByLocalID(V_ID(i - 1, j - 0, k - 0));
					verts[7] = mesh->NodeByLocalID(V_ID(i - 0, j - 0, k - 0));

					mesh->CreateCell(verts,nvf,numnodes,6);
				}
			}
		}
	}


	if( cut_grid )
	{
		for (Mesh::iteratorCell it = mesh->BeginCell(); it != mesh->EndCell(); ++it)
		{
			Storage::real cnt[3];
			it->Centroid(cnt);
			if( cnt[0] > 0.25 && cnt[0] < 0.75 && cnt[1] > 0.25 && cnt[1] < 0.75 )
				it->Delete();
			
		}

		for (Mesh::iteratorElement it = mesh->BeginElement(FACE|EDGE|NODE); it != mesh->EndElement(); ++it)
			if( it->nbAdjElements(CELL) == 0 ) it->Delete();
	}

	if( alpha > 0.0 ) //skewness
	{
		const double eps = 1.0e-6;
		for( Mesh::iteratorNode it = mesh->BeginNode(); it != mesh->EndNode(); ++it)
		{
			Storage::real_array coords = it->Coords();
			Storage::real h = 1.0e20;

			ElementArray<Cell> it_cells = it->getCells();

			for(int k = 0; k < it_cells.size(); ++k)
			{
				Storage::real maxmin[6];
				maxmin[0] = -1e20;
				maxmin[1] = 1e20;
				maxmin[2] = -1e20;
				maxmin[3] = 1e20;
				ElementArray<Node> nodes = it_cells[k].getNodes();
				for (ElementArray<Node>::iterator it = nodes.begin(); it != nodes.end(); it++)
				{

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
//.........这里部分代码省略.........

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

	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)
		if( !it->GetMarker(slice) && it->getBeg()->GetMarker(slice) && it->getEnd()->GetMarker(slice) )
		{
			it->SetMarker(slice);
			nmark++;
		}

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

	if( !Element::CheckConnectivity(&m) )
		std::cout << "Connectivity is broken" << std::endl;
		
	
	nslice = 0;
	MarkerType visited = m.CreateMarker();
	for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it)
	{
		ElementArray<Edge> edges = it->getEdges(slice);
		if( edges.size() >= 3 ) //these should form a triangle
		{
			//order edges
			ElementArray<Edge> order_edges(&m);
			order_edges.push_back(edges[0]);
			order_edges.SetMarker(visited);
			while(order_edges.size() != edges.size() )
			{
				for(int k = 0; k < edges.size(); ++k) if( !edges[k]->GetMarker(visited) )
				{
					if( edges[k]->getBeg() == order_edges.back()->getBeg() || edges[k]->getBeg() == order_edges.back()->getEnd() ||
						edges[k]->getEnd() == order_edges.back()->getBeg() || edges[k]->getEnd() == order_edges.back()->getEnd() )
					{
						order_edges.push_back(edges[k]);
						order_edges.back().SetMarker(visited);
					}
				}
			}
			edges.RemMarker(visited);
			Face f = m.CreateFace(order_edges).first;
			f.Bulk(sliced) = 1;
			Cell::SplitCell(it->self(),ElementArray<Face>(&m,1,f.GetHandle()),0);
			nslice++;
		}
	}
	m.ReleaseMarker(visited);

	std::cout << "sliced cells: " << nslice << std::endl;
	
	if( !Element::CheckConnectivity(&m) )
		std::cout << "Connectivity is broken" << std::endl;

	Tag material = m.CreateTag("MATERIAL",DATA_INTEGER,CELL,NONE,1);

	for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it)
	{
		double cnt[3];
		it->Centroid(cnt);
		double v = cnt[0]*nx+cnt[1]*ny+cnt[2]*nz-d;
		if( v < 0.0 )
			it->Integer(material) = 0;
		else
			it->Integer(material) = 1;
	}
	

	m.ReleaseMarker(slice,NODE|EDGE);

	m.Save(grid_out);
	return 0;
}