50 mMeshFilename(
"mesh"),
51 mDomainWidth(numNodesAlongWidth*scaleFactor),
52 mNumCellWidth(numNodesAlongWidth),
53 mNumCellLength(numNodesAlongLength)
59 assert(scaleFactor > 0.0);
62 std::stringstream pid;
64 OutputFileHandler output_file_handler(
"2D_temporary_honeycomb_mesh_" + pid.str());
70 double vertical_spacing = (sqrt(3.0)/2)*horizontal_spacing;
76 num_nodes_along_width += 2*ghosts;
77 num_nodes_along_length += 2*ghosts;
79 unsigned num_nodes = num_nodes_along_width*num_nodes_along_length;
80 unsigned num_elem_along_width = num_nodes_along_width-1;
81 unsigned num_elem_along_length = num_nodes_along_length-1;
82 unsigned num_elem = 2*num_elem_along_width*num_elem_along_length;
83 unsigned num_edges = 3*num_elem_along_width*num_elem_along_length + num_elem_along_width + num_elem_along_length;
85 double x0 = -horizontal_spacing*ghosts;
86 double y0 = -vertical_spacing*ghosts;
88 mBottom = -vertical_spacing*ghosts;
89 mTop =
mBottom + vertical_spacing*(num_nodes_along_length-1);
93 (*p_node_file) << std::scientific;
94 (*p_node_file) << std::setprecision(20);
95 (*p_node_file) << num_nodes <<
"\t2\t0\t1" << std::endl;
98 for (
unsigned i=0; i<num_nodes_along_length; i++)
100 for (
unsigned j=0; j<num_nodes_along_width; j++)
110 unsigned boundary = 0;
111 if (i==0 || i==num_nodes_along_length-1 || j==0 || j==num_nodes_along_width-1)
116 const double adjustment = i % 2 != 0 ? 0.5 : 0.0;
117 double x = x0 + horizontal_spacing*((
double)j + adjustment);
118 double y = y0 + vertical_spacing*(
double)i;
121 if ((y<0.0) && (y>-1e-12))
129 (*p_node_file) << node++ <<
"\t" << x <<
"\t" << y <<
"\t" << boundary << std::endl;
132 p_node_file->close();
136 (*p_elem_file) << std::scientific;
139 (*p_node_file) << std::scientific;
141 (*p_elem_file) << num_elem <<
"\t3\t0" << std::endl;
142 (*p_edge_file) << num_edges <<
"\t1" << std::endl;
146 for (
unsigned i=0; i<num_elem_along_length; i++)
148 for (
unsigned j=0; j < num_elem_along_width; j++)
150 unsigned node0 = i*num_nodes_along_width + j;
151 unsigned node1 = i*num_nodes_along_width + j+1;
152 unsigned node2 = (i+1)*num_nodes_along_width + j;
159 (*p_elem_file) << elem++ <<
"\t" << node0 <<
"\t" << node1 <<
"\t" << node2 << std::endl;
161 unsigned horizontal_edge_is_boundary_edge = 0;
162 unsigned vertical_edge_is_boundary_edge = 0;
165 horizontal_edge_is_boundary_edge = 1;
169 vertical_edge_is_boundary_edge = 1;
172 (*p_edge_file) << edge++ <<
"\t" << node0 <<
"\t" << node1 <<
"\t" << horizontal_edge_is_boundary_edge << std::endl;
173 (*p_edge_file) << edge++ <<
"\t" << node1 <<
"\t" << node2 <<
"\t" << 0 << std::endl;
174 (*p_edge_file) << edge++ <<
"\t" << node2 <<
"\t" << node0 <<
"\t" << vertical_edge_is_boundary_edge << std::endl;
176 node0 = i*num_nodes_along_width + j + 1;
182 node1 = (i+1)*num_nodes_along_width + j+1;
183 node2 = (i+1)*num_nodes_along_width + j;
185 (*p_elem_file) << elem++ <<
"\t" << node0 <<
"\t" << node1 <<
"\t" << node2 << std::endl;
189 for (
unsigned i=0; i<num_elem_along_length; i++)
191 unsigned node0, node1;
195 node0 = (i+1)*num_nodes_along_width - 1;
196 node1 = (i+2)*num_nodes_along_width - 1;
200 node0 = (i+1)*num_nodes_along_width;
201 node1 = (i)*num_nodes_along_width;
203 (*p_edge_file) << edge++ <<
"\t" << node0 <<
"\t" << node1 <<
"\t" << 1 << std::endl;
206 for (
unsigned j=0; j<num_elem_along_width; j++)
208 unsigned node0 = num_nodes_along_width*(num_nodes_along_length-1) + j;
209 unsigned node1 = num_nodes_along_width*(num_nodes_along_length-1) + j+1;
210 (*p_edge_file) << edge++ <<
"\t" << node1 <<
"\t" << node0 <<
"\t" << 1 << std::endl;
213 p_elem_file->close();
214 p_edge_file->close();
220 mpMesh = boost::make_shared<MutableMesh<2,2> >();
221 mpMesh->ConstructFromMeshReader(mesh_reader);
228 mpMesh->SetMeshHasChangedSinceLoading();
254 EXCEPTION(
"Cannot call GetCircularMesh on a HoneycombMesh with ghost nodes");
258 c_vector<double,2> centre = zero_vector<double>(2);
259 for (
unsigned i=0; i<
mpMesh->GetNumNodes(); i++)
261 centre +=
mpMesh->GetNode(i)->rGetLocation();
265 mpMesh->Translate(-centre[0], -centre[1]);
268 for (
unsigned i=0; i<
mpMesh->GetNumAllNodes(); i++)
270 if (norm_2(
mpMesh->GetNode(i)->rGetLocation()) >= radius)
272 mpMesh->DeleteNodePriorToReMesh(i);
277 c_vector<double,2>& r_location =
mpMesh->GetNode(i)->rGetModifiableLocation();
278 c_vector<double,2> shift;
280 double max_jiggle = radius*5e-6;
281 shift[0] = max_jiggle*(p_gen->
ranf()-0.5);
282 shift[1] = max_jiggle*(p_gen->
ranf()-0.5);