MultiSource/Benchmarks/DOE-ProxyApps-C++/miniFE/simple_mesh_description.hpp - third_party/llvm-test-suite - Git at Google


 #ifndef _simple_mesh_description_hpp_
 #define _simple_mesh_description_hpp_

 //@HEADER
 // ************************************************************************
 //
 // MiniFE: Simple Finite Element Assembly and Solve
 // Copyright (2006-2013) Sandia Corporation
 //
 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
 // license for use of this work by or on behalf of the U.S. Government.
 //
 // This library is free software; you can redistribute it and/or modify
 // it under the terms of the GNU Lesser General Public License as
 // published by the Free Software Foundation; either version 2.1 of the
 // License, or (at your option) any later version.
 //
 // This library is distributed in the hope that it will be useful, but
 // WITHOUT ANY WARRANTY; without even the implied warranty of
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Lesser General Public License for more details.
 //
 // You should have received a copy of the GNU Lesser General Public
 // License along with this library; if not, write to the Free Software
 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 // USA
 //
 // ************************************************************************
 //@HEADER

 #include <utils.hpp>
 #include <set>
 #include <map>

 namespace miniFE {

 template<typename GlobalOrdinal>
 class simple_mesh_description {
 public:
   simple_mesh_description(const Box& global_box_in, const Box& local_box_in)
   {
    Box local_node_box;
     for(int i=0; i<3; ++i) {
       global_box[i][0] = global_box_in[i][0];
       global_box[i][1] = global_box_in[i][1];
       local_box[i][0] = local_box_in[i][0];
       local_box[i][1] = local_box_in[i][1];
       local_node_box[i][0] = local_box_in[i][0];
       local_node_box[i][1] = local_box_in[i][1];
       //num-owned-nodes == num-elems+1 in this dimension if the elem box is not empty
       //and we are at the high end of the global range in that dimension:
       if (local_box_in[i][1] > local_box_in[i][0] && local_box_in[i][1] == global_box[i][1]) local_node_box[i][1] += 1;
     }

     int max_node_x = global_box[0][1]+1;
     int max_node_y = global_box[1][1]+1;
     int max_node_z = global_box[2][1]+1;
     create_map_id_to_row(max_node_x, max_node_y, max_node_z, local_node_box,
                          map_ids_to_rows);

     //As described in analytic_soln.hpp,
     //we will impose a 0 boundary-condition on faces x=0, y=0, z=0, y=1, z=1
     //we will impose a 1 boundary-condition on face x=1

 #ifdef MINIFE_DEBUG
 std::cout<<std::endl;
 #endif
     const int X=0;
     const int Y=1;
     const int Z=2;

     const int x1 = max_node_x - 1;
     const int y1 = max_node_y - 1;
     const int z1 = max_node_z - 1;

     //if we're on the x=0 face:
     if (global_box[X][0] == local_box[X][0]) {
       int miny = local_node_box[Y][0];
       int minz = local_node_box[Z][0];
       int maxy = local_node_box[Y][1];
       int maxz = local_node_box[Z][1];
       //expand y and z dimensions to include ghost layer
       if (local_node_box[Y][0] > 0) --miny;
       if (local_node_box[Z][0] > 0) --minz;
       if (local_node_box[Y][1] < max_node_y) ++maxy;
       if (local_node_box[Z][1] < max_node_z) ++maxz;

       for(int iz=minz; iz<maxz; ++iz) {
         for(int iy=miny; iy<maxy; ++iy) {
           GlobalOrdinal nodeID = get_id<GlobalOrdinal>(max_node_x, max_node_y, max_node_z,
              0, iy, iz);
 #ifdef MINIFE_DEBUG
 std::cout<<"x=0 BC, node "<<nodeID<<", (0,"<<iy<<","<<iz<<")"<<std::endl;
 #endif
           bc_rows_0.insert(map_id_to_row(nodeID));
         }
       }
     }

     //if we're on the y=0 face:
     if (global_box[Y][0] == local_box[Y][0]) {
       int minx = local_node_box[X][0];
       int minz = local_node_box[Z][0];
       int maxx = local_node_box[X][1];
       int maxz = local_node_box[Z][1];
       //expand x and z dimensions to include ghost layer
       if (local_node_box[X][0] > 0) --minx;
       if (local_node_box[Z][0] > 0) --minz;
       if (local_node_box[X][1] < max_node_x) ++maxx;
       if (local_node_box[Z][1] < max_node_z) ++maxz;

       for(int iz=minz; iz<maxz; ++iz) {
         for(int ix=minx; ix<maxx; ++ix) {
           GlobalOrdinal nodeID = get_id<GlobalOrdinal>(max_node_x, max_node_y, max_node_z,
              ix, 0, iz);
 #ifdef MINIFE_DEBUG
 std::cout<<"y=0 BC, node "<<nodeID<<", ("<<ix<<",0,"<<iz<<")"<<std::endl;
 #endif
           GlobalOrdinal row = map_id_to_row(nodeID);
           if (row < 0) {
             std::cout<<"on the y==0 face (ix="<<ix<<", iz="<<iz<<"), ERROR: found negative row ("<<row<<") for nodeID="<<nodeID<<std::endl;
           }
           bc_rows_0.insert(row);
         }
       }
     }

     //if we're on the z=0 face:
     if (global_box[Z][0] == local_box[Z][0]) {
       int minx = local_node_box[X][0];
       int miny = local_node_box[Y][0];
       int maxx = local_node_box[X][1];
       int maxy = local_node_box[Y][1];
       //expand x and y dimensions to include ghost layer
       if (local_node_box[X][0] > 0) --minx;
       if (local_node_box[Y][0] > 0) --miny;
       if (local_node_box[X][1] < max_node_x) ++maxx;
       if (local_node_box[Y][1] < max_node_y) ++maxy;

       for(int iy=miny; iy<maxy; ++iy) {
         for(int ix=minx; ix<maxx; ++ix) {
           GlobalOrdinal nodeID = get_id<GlobalOrdinal>(max_node_x, max_node_y, max_node_z,
              ix, iy, 0);
 #ifdef MINIFE_DEBUG
 std::cout<<"z=0 BC, node "<<nodeID<<", ("<<ix<<","<<iy<<",0)"<<std::endl;
 #endif
           bc_rows_0.insert(map_id_to_row(nodeID));
         }
       }
     }

     //if we're on the x=1 face:
     if (global_box[X][1] == local_box[X][1]) {
       int minz = local_node_box[Z][0];
       int miny = local_node_box[Y][0];
       int maxz = local_node_box[Z][1];
       int maxy = local_node_box[Y][1];
       //expand z and y dimensions to include ghost layer
       if (local_node_box[Z][0] > 0) --minz;
       if (local_node_box[Y][0] > 0) --miny;
       if (local_node_box[Z][1] < max_node_z) ++maxz;
       if (local_node_box[Y][1] < max_node_y) ++maxy;

       for(int iy=miny; iy<maxy; ++iy) {
         for(int iz=minz; iz<maxz; ++iz) {
           GlobalOrdinal nodeID = get_id<GlobalOrdinal>(max_node_x, max_node_y, max_node_z,
              x1, iy, iz);
           GlobalOrdinal row = map_id_to_row(nodeID);
 #ifdef MINIFE_DEBUG
 std::cout<<"x=1 BC, node "<<nodeID<<", row "<<row<<", ("<<x1<<","<<iy<<","<<iz<<")"<<std::endl;
 #endif
           bc_rows_1.insert(row);
         }
       }
     }

     //if we're on the y=1 face:
     if (global_box[Y][1] == local_box[Y][1]) {
       int minz = local_node_box[Z][0];
       int minx = local_node_box[X][0];
       int maxz = local_node_box[Z][1];
       int maxx = local_node_box[X][1];
       //expand z and x dimensions to include ghost layer
       if (local_node_box[Z][0] > 0) --minz;
       if (local_node_box[X][0] > 0) --minx;
       if (local_node_box[Z][1] < max_node_z) ++maxz;
       if (local_node_box[X][1] < max_node_x) ++maxx;

       for(int ix=minx; ix<maxx; ++ix) {
         for(int iz=minz; iz<maxz; ++iz) {
           GlobalOrdinal nodeID = get_id<GlobalOrdinal>(max_node_x, max_node_y, max_node_z,
              ix, y1, iz);
 #ifdef MINIFE_DEBUG
 std::cout<<"y=1 BC, node "<<nodeID<<", ("<<ix<<","<<y1<<","<<iz<<")"<<std::endl;
 #endif
           bc_rows_0.insert(map_id_to_row(nodeID));
         }
       }
     }

     //if we're on the z=1 face:
     if (global_box[Z][1] == local_box[Z][1]) {
       int miny = local_node_box[Y][0];
       int minx = local_node_box[X][0];
       int maxy = local_node_box[Y][1];
       int maxx = local_node_box[X][1];
       //expand x and y dimensions to include ghost layer
       if (local_node_box[Y][0] > 0) --miny;
       if (local_node_box[X][0] > 0) --minx;
       if (local_node_box[Y][1] < max_node_y) ++maxy;
       if (local_node_box[X][1] < max_node_x) ++maxx;

       for(int ix=minx; ix<maxx; ++ix) {
         for(int iy=miny; iy<maxy; ++iy) {
           GlobalOrdinal nodeID = get_id<GlobalOrdinal>(max_node_x, max_node_y, max_node_z,
              ix, iy, z1);
 #ifdef MINIFE_DEBUG
 std::cout<<"z=1 BC, node "<<nodeID<<", ("<<ix<<","<<iy<<","<<z1<<")"<<std::endl;
 #endif
           bc_rows_0.insert(map_id_to_row(nodeID));
         }
       }
     }

   }

   GlobalOrdinal map_id_to_row(const GlobalOrdinal& id) const
   {
     return find_row_for_id(id, map_ids_to_rows);
   }

   GlobalOrdinal max_row_in_map() const {
     if (map_ids_to_rows.empty()) return 0;
     typename std::map<GlobalOrdinal,GlobalOrdinal>::const_iterator mend = map_ids_to_rows.end();
     --mend;
     return mend->second;
   }
   std::set<GlobalOrdinal> bc_rows_0;
   std::set<GlobalOrdinal> bc_rows_1;
   std::map<GlobalOrdinal,GlobalOrdinal> map_ids_to_rows;
   Box global_box;
   Box local_box;
 };//class simple_mesh_description

 }//namespace miniFE

 #endif

	#ifndef _simple_mesh_description_hpp_
	#define _simple_mesh_description_hpp_

	//@HEADER
	// ************************************************************************
	//
	// MiniFE: Simple Finite Element Assembly and Solve
	// Copyright (2006-2013) Sandia Corporation
	//
	// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
	// license for use of this work by or on behalf of the U.S. Government.
	//
	// This library is free software; you can redistribute it and/or modify
	// it under the terms of the GNU Lesser General Public License as
	// published by the Free Software Foundation; either version 2.1 of the
	// License, or (at your option) any later version.
	//
	// This library is distributed in the hope that it will be useful, but
	// WITHOUT ANY WARRANTY; without even the implied warranty of
	// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
	// Lesser General Public License for more details.
	//
	// You should have received a copy of the GNU Lesser General Public
	// License along with this library; if not, write to the Free Software
	// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
	// USA
	//
	// ************************************************************************
	//@HEADER

	#include <utils.hpp>
	#include <set>
	#include <map>

	namespace miniFE {

	template<typename GlobalOrdinal>
	class simple_mesh_description {
	public:
	simple_mesh_description(const Box& global_box_in, const Box& local_box_in)
	{
	Box local_node_box;
	for(int i=0; i<3; ++i) {
	global_box[i][0] = global_box_in[i][0];
	global_box[i][1] = global_box_in[i][1];
	local_box[i][0] = local_box_in[i][0];
	local_box[i][1] = local_box_in[i][1];
	local_node_box[i][0] = local_box_in[i][0];
	local_node_box[i][1] = local_box_in[i][1];
	//num-owned-nodes == num-elems+1 in this dimension if the elem box is not empty
	//and we are at the high end of the global range in that dimension:
	if (local_box_in[i][1] > local_box_in[i][0] && local_box_in[i][1] == global_box[i][1]) local_node_box[i][1] += 1;
	}

	int max_node_x = global_box[0][1]+1;
	int max_node_y = global_box[1][1]+1;
	int max_node_z = global_box[2][1]+1;
	create_map_id_to_row(max_node_x, max_node_y, max_node_z, local_node_box,
	map_ids_to_rows);

	//As described in analytic_soln.hpp,
	//we will impose a 0 boundary-condition on faces x=0, y=0, z=0, y=1, z=1
	//we will impose a 1 boundary-condition on face x=1

	#ifdef MINIFE_DEBUG
	std::cout<<std::endl;
	#endif
	const int X=0;
	const int Y=1;
	const int Z=2;

	const int x1 = max_node_x - 1;
	const int y1 = max_node_y - 1;
	const int z1 = max_node_z - 1;

	//if we're on the x=0 face:
	if (global_box[X][0] == local_box[X][0]) {
	int miny = local_node_box[Y][0];
	int minz = local_node_box[Z][0];
	int maxy = local_node_box[Y][1];
	int maxz = local_node_box[Z][1];
	//expand y and z dimensions to include ghost layer
	if (local_node_box[Y][0] > 0) --miny;
	if (local_node_box[Z][0] > 0) --minz;
	if (local_node_box[Y][1] < max_node_y) ++maxy;
	if (local_node_box[Z][1] < max_node_z) ++maxz;

	for(int iz=minz; iz<maxz; ++iz) {
	for(int iy=miny; iy<maxy; ++iy) {
	GlobalOrdinal nodeID = get_id<GlobalOrdinal>(max_node_x, max_node_y, max_node_z,
	0, iy, iz);
	#ifdef MINIFE_DEBUG
	std::cout<<"x=0 BC, node "<<nodeID<<", (0,"<<iy<<","<<iz<<")"<<std::endl;
	#endif
	bc_rows_0.insert(map_id_to_row(nodeID));
	}
	}
	}

	//if we're on the y=0 face:
	if (global_box[Y][0] == local_box[Y][0]) {
	int minx = local_node_box[X][0];
	int minz = local_node_box[Z][0];
	int maxx = local_node_box[X][1];
	int maxz = local_node_box[Z][1];
	//expand x and z dimensions to include ghost layer
	if (local_node_box[X][0] > 0) --minx;
	if (local_node_box[Z][0] > 0) --minz;
	if (local_node_box[X][1] < max_node_x) ++maxx;
	if (local_node_box[Z][1] < max_node_z) ++maxz;

	for(int iz=minz; iz<maxz; ++iz) {
	for(int ix=minx; ix<maxx; ++ix) {
	GlobalOrdinal nodeID = get_id<GlobalOrdinal>(max_node_x, max_node_y, max_node_z,
	ix, 0, iz);
	#ifdef MINIFE_DEBUG
	std::cout<<"y=0 BC, node "<<nodeID<<", ("<<ix<<",0,"<<iz<<")"<<std::endl;
	#endif
	GlobalOrdinal row = map_id_to_row(nodeID);
	if (row < 0) {
	std::cout<<"on the y==0 face (ix="<<ix<<", iz="<<iz<<"), ERROR: found negative row ("<<row<<") for nodeID="<<nodeID<<std::endl;
	}
	bc_rows_0.insert(row);
	}
	}
	}

	//if we're on the z=0 face:
	if (global_box[Z][0] == local_box[Z][0]) {
	int minx = local_node_box[X][0];
	int miny = local_node_box[Y][0];
	int maxx = local_node_box[X][1];
	int maxy = local_node_box[Y][1];
	//expand x and y dimensions to include ghost layer
	if (local_node_box[X][0] > 0) --minx;
	if (local_node_box[Y][0] > 0) --miny;
	if (local_node_box[X][1] < max_node_x) ++maxx;
	if (local_node_box[Y][1] < max_node_y) ++maxy;

	for(int iy=miny; iy<maxy; ++iy) {
	for(int ix=minx; ix<maxx; ++ix) {
	GlobalOrdinal nodeID = get_id<GlobalOrdinal>(max_node_x, max_node_y, max_node_z,
	ix, iy, 0);
	#ifdef MINIFE_DEBUG
	std::cout<<"z=0 BC, node "<<nodeID<<", ("<<ix<<","<<iy<<",0)"<<std::endl;
	#endif
	bc_rows_0.insert(map_id_to_row(nodeID));
	}
	}
	}

	//if we're on the x=1 face:
	if (global_box[X][1] == local_box[X][1]) {
	int minz = local_node_box[Z][0];
	int miny = local_node_box[Y][0];
	int maxz = local_node_box[Z][1];
	int maxy = local_node_box[Y][1];
	//expand z and y dimensions to include ghost layer
	if (local_node_box[Z][0] > 0) --minz;
	if (local_node_box[Y][0] > 0) --miny;
	if (local_node_box[Z][1] < max_node_z) ++maxz;
	if (local_node_box[Y][1] < max_node_y) ++maxy;

	for(int iy=miny; iy<maxy; ++iy) {
	for(int iz=minz; iz<maxz; ++iz) {
	GlobalOrdinal nodeID = get_id<GlobalOrdinal>(max_node_x, max_node_y, max_node_z,
	x1, iy, iz);
	GlobalOrdinal row = map_id_to_row(nodeID);
	#ifdef MINIFE_DEBUG
	std::cout<<"x=1 BC, node "<<nodeID<<", row "<<row<<", ("<<x1<<","<<iy<<","<<iz<<")"<<std::endl;
	#endif
	bc_rows_1.insert(row);
	}
	}
	}

	//if we're on the y=1 face:
	if (global_box[Y][1] == local_box[Y][1]) {
	int minz = local_node_box[Z][0];
	int minx = local_node_box[X][0];
	int maxz = local_node_box[Z][1];
	int maxx = local_node_box[X][1];
	//expand z and x dimensions to include ghost layer
	if (local_node_box[Z][0] > 0) --minz;
	if (local_node_box[X][0] > 0) --minx;
	if (local_node_box[Z][1] < max_node_z) ++maxz;
	if (local_node_box[X][1] < max_node_x) ++maxx;

	for(int ix=minx; ix<maxx; ++ix) {
	for(int iz=minz; iz<maxz; ++iz) {
	GlobalOrdinal nodeID = get_id<GlobalOrdinal>(max_node_x, max_node_y, max_node_z,
	ix, y1, iz);
	#ifdef MINIFE_DEBUG
	std::cout<<"y=1 BC, node "<<nodeID<<", ("<<ix<<","<<y1<<","<<iz<<")"<<std::endl;
	#endif
	bc_rows_0.insert(map_id_to_row(nodeID));
	}
	}
	}

	//if we're on the z=1 face:
	if (global_box[Z][1] == local_box[Z][1]) {
	int miny = local_node_box[Y][0];
	int minx = local_node_box[X][0];
	int maxy = local_node_box[Y][1];
	int maxx = local_node_box[X][1];
	//expand x and y dimensions to include ghost layer
	if (local_node_box[Y][0] > 0) --miny;
	if (local_node_box[X][0] > 0) --minx;
	if (local_node_box[Y][1] < max_node_y) ++maxy;
	if (local_node_box[X][1] < max_node_x) ++maxx;

	for(int ix=minx; ix<maxx; ++ix) {
	for(int iy=miny; iy<maxy; ++iy) {
	GlobalOrdinal nodeID = get_id<GlobalOrdinal>(max_node_x, max_node_y, max_node_z,
	ix, iy, z1);
	#ifdef MINIFE_DEBUG
	std::cout<<"z=1 BC, node "<<nodeID<<", ("<<ix<<","<<iy<<","<<z1<<")"<<std::endl;
	#endif
	bc_rows_0.insert(map_id_to_row(nodeID));
	}
	}
	}

	}

	GlobalOrdinal map_id_to_row(const GlobalOrdinal& id) const
	{
	return find_row_for_id(id, map_ids_to_rows);
	}

	GlobalOrdinal max_row_in_map() const {
	if (map_ids_to_rows.empty()) return 0;
	typename std::map<GlobalOrdinal,GlobalOrdinal>::const_iterator mend = map_ids_to_rows.end();
	--mend;
	return mend->second;
	}
	std::set<GlobalOrdinal> bc_rows_0;
	std::set<GlobalOrdinal> bc_rows_1;
	std::map<GlobalOrdinal,GlobalOrdinal> map_ids_to_rows;
	Box global_box;
	Box local_box;
	};//class simple_mesh_description

	}//namespace miniFE

	#endif