
// Now we can either modify, add or delete existing surfaces and regions

// Let us first override the characteristic length parameter
lc = 10 ;

// (meshx) X (meshy) 
meshx = 5 ;
meshy = 5 ;

x_size = 15;
y_size = 15;

// Now let us include the 2D lattice assembly geometry routine ;
// And create the lattice points, lines and surfaces

Printf ('Characteristic length = %g', lc) ;
Printf ('Mesh X = %g, Y = %g', meshx, meshy) ;
Printf ('Ass X = %g, Y = %g', x_size, y_size) ;

// evenly distribute the points based on num_x, num_y and sizes
pointnum = 0 ;

Printf ('Generating points.') ;

For nx In {1:meshx+1}

	For ny In {1:meshy+1}
	
		pointnum = pointnum + 1 ;
		Point(pointnum) = {x_size*(nx-1), x_size*(ny-1), 0, lc};

	EndFor
	
EndFor

Printf ('Finished generating points.') ;

Printf ('Build lines based on points.') ;
// connect the points with lines to create a square grid
linenum = pointnum  ;
stub1 = linenum ;
// connect points in x direction first
For nx In {1:meshx+1}

	For ny In {1:meshy}
	
		linenum = linenum + 1 ;
		Line(linenum) = {(nx-1)*(meshy+1)+ny,(nx-1)*(meshy+1)+ny+1} ;

	EndFor
	
EndFor
stub2 = linenum ;
// now connect points in y direction
For ny In {1:meshy+1}

	For nx In {1:meshx}

		linenum = linenum + 1 ;
		Line(linenum) = {(nx-1)*(meshy+1)+ny,(nx)*(meshy+1)+ny} ;

	EndFor
	
EndFor
Printf ('Finished generating lines that connects the points.') ;

// let us now associate the lines together to create surfaces
surfnum = linenum ;
surfstride = 2 ;

isurf = 1 ;
For nx In {1:meshx}

	For ny In {1:meshy}

		surfnum = surfnum + 1 ;
		// Printf("Surface %g  = [%g,%g,%g,%g]",surfnum,-(stub1+(nx-1)*(meshy+1)+ny)+nx-1,stub2+(ny-1)*(meshx+1)+nx-ny+1,stub1+(nx-1)*(meshy+1)+ny-nx+1+meshy,-(stub2+(nx)*(meshy+1)+ny-1+meshx)) ;

		Line Loop(surfnum) = {-(stub1+(nx-1)*(meshy+1)+ny)+nx-1,stub2+(ny-1)*(meshx+1)+nx-ny+1,stub1+(nx-1)*(meshy+1)+ny-nx+1+meshy,-(stub2+(ny-1)*(meshx+1)+nx-ny+1+meshx)} ;
		Plane Surface(surfnum+1) = {surfnum} ;
		
		Physical Surface(Sprintf("Region%02g", isurf)) = { surfnum+1 } ;
		//Physical Surface(Sprintf("Region%05g", isurf+100)) = { surfnum } ;
		isurf = isurf + 1 ;
		
		surfnum = surfnum + 1 ;

	EndFor
	
EndFor

Printf ('Finished generating surfaces based on lines.') ;
