Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

source: issm/trunk/src/c/classes/AmrNeopz.cpp@ 25836

Last change on this file since 25836 was 25836, checked in by Mathieu Morlighem, 4 years ago
merged trunk-jpl and trunk for revision 25834
File size: 34.8 KB

Rev	Line
[21484]	1	/*!\file AdaptiveMeshrefinement.cpp
	2	* \brief: implementation of the adaptive mesh refinement tool based on NeoPZ library: github.com/labmec/neopz
	3	*/
	4
[21528]	5	#ifdef HAVE_CONFIG_H
[21484]	6	#include <config.h>
	7	#else
	8	#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
	9	#endif
	10
[25454]	11	#include "./AmrNeopz.h"
[21919]	12
[24496]	13	/Includes/
	14	/{{{/
	15	/Common includes/
	16	#include <iostream>
	17	#include <fstream>
	18	#include <string>
	19	#include <climits>
	20	#include <cfloat>
	21
	22	/NeoPZ includes/
	23	#include <pz_config.h>
	24	#include <pzreal.h>
	25	#include <pzvec.h>
	26	#include <pzeltype.h>
	27
	28	#include <TPZRefPatternTools.h>
	29	#include <TPZRefPatternDataBase.h>
	30	#include <TPZRefPattern.h>
	31
	32	#include <tpzchangeel.h>
	33	#include <TPZGeoElement.h>
	34	#include <pzreftriangle.h>
	35	#include <pzgeotriangle.h>
	36	#include <tpzgeoelrefpattern.h>
	37	#include <pzgraphmesh.h>
	38	#include <TPZVTKGeoMesh.h>
	39	/}}}/
	40
[21806]	41	using namespace pzgeom;
[21484]	42
	43	/Constructor, copy, clean up and destructor/
[25454]	44	AmrNeopz::AmrNeopz(){/{{{/
[23501]	45
	46	/Set pointers to NULL/
	47	this->fathermesh = NULL;
	48	this->previousmesh = NULL;
	49	this->refinement_type = -1;
	50	this->level_max = -1;
	51	this->gradation = -1;
	52	this->lag = -1;
	53	this->groundingline_distance = -1;
	54	this->icefront_distance = -1;
	55	this->thicknesserror_threshold = -1;
	56	this->deviatoricerror_threshold = -1;
	57	this->deviatoricerror_maximum = -1;
	58	this->thicknesserror_maximum = -1;
	59	this->sid2index.clear();
	60	this->index2sid.clear();
	61	this->specialelementsindex.clear();
	62	this->x = NULL;
	63	this->y = NULL;
	64	this->elementslist = NULL;
	65	this->numberofvertices = -1;
	66	this->numberofelements = -1;
[21484]	67	}
	68	/}}}/
[25454]	69	AmrNeopz::AmrNeopz(const AmrNeopz &cp){/{{{/
[21806]	70	this->Initialize();
[21562]	71	this->operator =(cp);
[21484]	72	}
	73	/}}}/
[25454]	74	AmrNeopz & AmrNeopz::operator =(const AmrNeopz &cp){/{{{/
[21484]	75
[21562]	76	/Clean all attributes/
	77	this->CleanUp();
	78	/Copy all data/
[22100]	79	this->fathermesh = new TPZGeoMesh(*cp.fathermesh);
	80	this->previousmesh = new TPZGeoMesh(*cp.previousmesh);
[22241]	81	this->refinement_type = cp.refinement_type;
[22100]	82	this->level_max = cp.level_max;
	83	this->gradation = cp.gradation;
	84	this->lag = cp.lag;
	85	this->groundingline_distance = cp.groundingline_distance;
	86	this->icefront_distance = cp.icefront_distance;
	87	this->thicknesserror_threshold = cp.thicknesserror_threshold;
	88	this->deviatoricerror_threshold = cp.deviatoricerror_threshold;
[22276]	89	this->deviatoricerror_maximum = cp.deviatoricerror_maximum;
	90	this->thicknesserror_maximum = cp.thicknesserror_maximum;
[21806]	91	this->sid2index.clear();
	92	this->sid2index.resize(cp.sid2index.size());
	93	for(int i=0;i<cp.sid2index.size();i++) this->sid2index[i]=cp.sid2index[i];
[21919]	94	this->index2sid.clear();
	95	this->index2sid.resize(cp.index2sid.size());
	96	for(int i=0;i<cp.index2sid.size();i++) this->index2sid[i]=cp.index2sid[i];
	97	this->specialelementsindex.clear();
	98	this->specialelementsindex.resize(cp.specialelementsindex.size());
	99	for(int i=0;i<cp.specialelementsindex.size();i++) this->specialelementsindex[i]=cp.specialelementsindex[i];
	100
[21562]	101	return *this;
[21484]	102	}
	103	/}}}/
[25454]	104	AmrNeopz::~AmrNeopz(){/{{{/
[23501]	105	int writemesh = 0;//only to restart
[23471]	106	if(writemesh) this->WriteMesh();
[21562]	107	this->CleanUp();
	108	gRefDBase.clear();
[21484]	109	}
	110	/}}}/
[25454]	111	void AmrNeopz::CleanUp(){/{{{/
[21484]	112
[21806]	113	/Verify and delete all data/
[21562]	114	if(this->fathermesh) delete this->fathermesh;
[22100]	115	if(this->previousmesh) delete this->previousmesh;
[24499]	116	if(this->x) xDelete<IssmDouble>(this->x);
	117	if(this->y) xDelete<IssmDouble>(this->y);
	118	if(this->elementslist) xDelete<int>(this->elementslist);
[22241]	119	this->refinement_type = -1;
[22100]	120	this->level_max = -1;
	121	this->gradation = -1;
	122	this->lag = -1;
	123	this->groundingline_distance = -1;
	124	this->icefront_distance = -1;
	125	this->thicknesserror_threshold = -1;
	126	this->deviatoricerror_threshold = -1;
[22276]	127	this->deviatoricerror_maximum = -1;
	128	this->thicknesserror_maximum = -1;
[23501]	129	this->numberofvertices = -1;
	130	this->numberofelements = -1;
[21806]	131	this->sid2index.clear();
[21919]	132	this->index2sid.clear();
	133	this->specialelementsindex.clear();
[21484]	134	}
	135	/}}}/
[21562]	136
[23501]	137	/Mesh refinement methods/
[25454]	138	void AmrNeopz::SetMesh(int elementslist_in,IssmDouble x_in,IssmDouble** y_in,int* numberofvertices_in,int* numberofelements_in){/{{{/
[23501]	139
	140	/Delete previous mesh and keep the entire mesh/
	141	if(this->elementslist) xDelete<int>(this->elementslist);
	142	if(this->x) xDelete<IssmDouble>(this->x);
	143	if(this->y) xDelete<IssmDouble>(this->y);
[21484]	144
[23501]	145	this->elementslist = *elementslist_in;
	146	this->x = *x_in;
	147	this->y = *y_in;
	148	this->numberofvertices = *numberofvertices_in;
	149	this->numberofelements = *numberofelements_in;
	150	}/}}}/
[25454]	151	void AmrNeopz::GetMesh(int elementslist_out,IssmDouble x_out,IssmDouble** y_out,int* numberofvertices_out,int* numberofelements_out){/{{{/
[23501]	152
	153	/Get the entire mesh/
	154	*elementslist_out = this->elementslist;
	155	*x_out = this->x;
	156	*y_out = this->y;
	157	*numberofvertices_out= this->numberofvertices;
	158	*numberofelements_out= this->numberofelements;
	159	}/}}}/
[25454]	160	void AmrNeopz::ExecuteRefinement(double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror,int pdatalist,double pxylist,int** pelementslist){/{{{/
[23066]	161
[22100]	162	/IMPORTANT! pelementslist are in Matlab indexing/
[21528]	163	/NEOPZ works only in C indexing/
[22100]	164	if(!this->fathermesh \|\| !this->previousmesh) _error_("Impossible to execute refinement: fathermesh or previousmesh is NULL!\n");
[22241]	165	if(this->refinement_type!=0 && this->refinement_type!=1) _error_("Impossible to execute refinement: refinement type is not defined!\n");
[21528]	166
[22241]	167	/Input verifications/
	168	if(this->deviatoricerror_threshold>0 && !deviatoricerror) _error_("deviatoricerror is NULL!\n");
	169	if(this->thicknesserror_threshold>0 && !thicknesserror) _error_("thicknesserror is NULL!\n");
	170	if(this->groundingline_distance>0 && !gl_distance) _error_("gl_distance is NULL!\n");
	171	if(this->icefront_distance>0 && !if_distance) _error_("if_distance is NULL!\n");
[22294]	172	/Attributes verifications/
	173	if(this->deviatoricerror_threshold>0 && this->deviatoricerror_groupthreshold<DBL_EPSILON) _error_("group threshold is too small!");
	174	if(this->thicknesserror_threshold>0 && this->thicknesserror_groupthreshold<DBL_EPSILON) _error_("group threshold is too small!");
[22241]	175
[22100]	176	/Intermediaries/
	177	bool verbose=VerboseSolution();
[23066]	178
[22100]	179	/Execute refinement/
[22241]	180	this->RefineMeshOneLevel(verbose,gl_distance,if_distance,deviatoricerror,thicknesserror);
[23066]	181
[21806]	182	/Get new geometric mesh in ISSM data structure/
[23065]	183	this->GetMesh(this->previousmesh,pdatalist,pxylist,pelementslist);
[21919]	184
[21806]	185	/Verify the new geometry/
[23065]	186	this->CheckMesh(pdatalist,pxylist,pelementslist);
[21484]	187	}
	188	/}}}/
[25454]	189	void AmrNeopz::RefineMeshOneLevel(bool &verbose,double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror){/{{{/
[23066]	190
[21806]	191	/Intermediaries/
[22241]	192	int nelem =-1;
	193	int side2D = 6;
	194	int sid =-1;
	195	int count =-1;
	196	int criteria =-1;
	197	int numberofcriteria =-1;
	198	int nconformelements = this->sid2index.size();
[22276]	199	double gl_distance_h =-1;
	200	double gl_distance_hmax = this->groundingline_distance;
	201	double if_distance_h =-1;
	202	double if_distance_hmax = this->icefront_distance;
[22241]	203	double gl_groupdistance =-1;
	204	double if_groupdistance =-1;
	205	double d_maxerror =-1;
	206	double t_maxerror =-1;
	207	double deviatoric_grouperror =-1;
	208	double thickness_grouperror =-1;
	209	TPZGeoMesh* gmesh = NULL;
	210	TPZVec<REAL> qsi(2,0.),cp(3,0.);
	211	TPZVec<TPZGeoEl*> sons;
	212	std::vector<int> index;
[21919]	213
[22241]	214	/Calculate the number of criteria{{{/
	215	numberofcriteria=0;
	216	if(this->deviatoricerror_threshold>0) numberofcriteria++;
	217	if(this->thicknesserror_threshold>0) numberofcriteria++;
	218	if(this->groundingline_distance>0) numberofcriteria++;
	219	if(this->icefront_distance>0) numberofcriteria++;
	220	/}}}/
[21919]	221
[22241]	222	/Calculate the maximum of the estimators, if requested{{{/
[22276]	223	if(this->deviatoricerror_threshold>0 && this->deviatoricerror_maximum<DBL_EPSILON){
	224	for(int i=0;i<nconformelements;i++) this->deviatoricerror_maximum=max(this->deviatoricerror_maximum,deviatoricerror[i]);
[21806]	225	}
[22276]	226	if(this->thicknesserror_threshold>0 && this->thicknesserror_maximum<DBL_EPSILON){
	227	for(int i=0;i<nconformelements;i++) this->thicknesserror_maximum=max(this->thicknesserror_maximum,thicknesserror[i]);
[22241]	228	}
	229	/}}}/
	230
	231	/First, verify if special elements have min distance or high errors{{{/
	232	gmesh=this->previousmesh;
[21919]	233	for(int i=0;i<this->specialelementsindex.size();i++){
[22241]	234	if(this->specialelementsindex[i]==-1) _error_("index is -1!\n");
	235	if(!gmesh->Element(this->specialelementsindex[i])) _error_("element is null!\n");
	236	if(!gmesh->Element(this->specialelementsindex[i])->Father()) _error_("father is null!\n");
	237	if(gmesh->Element(this->specialelementsindex[i])->HasSubElement()) _error_("special element has sub elements!\n");
	238	sons.clear();
	239	gmesh->Element(this->specialelementsindex[i])->Father()->GetHigherSubElements(sons);
	240	/Limits/
[22276]	241	gl_distance_h = gl_distance_hmax*std::pow(this->gradation,this->level_max-gmesh->Element(this->specialelementsindex[i])->Level());
	242	if_distance_h = if_distance_hmax*std::pow(this->gradation,this->level_max-gmesh->Element(this->specialelementsindex[i])->Level());
	243	d_maxerror = this->deviatoricerror_threshold*this->deviatoricerror_maximum;
	244	t_maxerror = this->thicknesserror_threshold*this->thicknesserror_maximum;
[22241]	245	/Calculate the distance and error of the group (sons)/
	246	gl_groupdistance=INFINITY;if_groupdistance=INFINITY;deviatoric_grouperror=0;thickness_grouperror=0;
	247	for(int s=0;s<sons.size();s++){
	248	sid=this->index2sid[sons[s]->Index()];
	249	if(sid<0) continue;
	250	if(this->groundingline_distance>0) gl_groupdistance=std::min(gl_groupdistance,gl_distance[sid]);
	251	if(this->icefront_distance>0) if_groupdistance=std::min(if_groupdistance,if_distance[sid]);
	252	if(this->deviatoricerror_threshold>0) deviatoric_grouperror+=deviatoricerror[sid];
	253	if(this->thicknesserror_threshold>0) thickness_grouperror+=thicknesserror[sid];
	254	}
	255	criteria=0;
[22276]	256	if(this->groundingline_distance>0 && gl_groupdistance<gl_distance_h+DBL_EPSILON) criteria++;
	257	if(this->icefront_distance>0 && if_groupdistance<if_distance_h+DBL_EPSILON) criteria++;
	258	if(this->deviatoricerror_threshold>0 && deviatoric_grouperror>d_maxerror-DBL_EPSILON) criteria++;
	259	if(this->thicknesserror_threshold>0 && thickness_grouperror>t_maxerror-DBL_EPSILON) criteria++;
[22241]	260	/Finally, it keeps the father index if it must be refine/
	261	if(criteria) index.push_back(gmesh->Element(this->specialelementsindex[i])->FatherIndex());
[21919]	262	}
[22241]	263	/}}}/
[23066]	264
[22241]	265	/Now, detele the special elements{{{/
	266	if(this->refinement_type==1) this->DeleteSpecialElements(verbose,gmesh);
	267	else this->specialelementsindex.clear();
	268	/}}}/
[21919]	269
[22241]	270	/Set the mesh and delete previousmesh if refinement type is 0{{{/
	271	if(this->refinement_type==0){
	272	delete this->previousmesh;
	273	gmesh=this->fathermesh;
	274	}
	275	/}}}/
[23066]	276
[22241]	277	/Unrefinement process: loop over level of refinements{{{/
	278	if(verbose) _printf_("\tunrefinement process...\n");
	279	if(verbose) _printf_("\ttotal: ");
	280	count=0;
[23066]	281
[22241]	282	nelem=gmesh->NElements();//must keep here
	283	for(int i=0;i<nelem;i++){
	284	if(!gmesh->Element(i)) continue;
	285	if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
	286	if(gmesh->Element(i)->HasSubElement()) continue;
	287	if(gmesh->Element(i)->Level()==0) continue;
	288	if(!gmesh->Element(i)->Father()) _error_("father is NULL!\n");
	289	/Limits with lag/
[22276]	290	gl_distance_h = this->laggl_distance_hmaxstd::pow(this->gradation,this->level_max-gmesh->Element(i)->Level());
	291	if_distance_h = this->lagif_distance_hmaxstd::pow(this->gradation,this->level_max-gmesh->Element(i)->Level());
[22294]	292	d_maxerror = this->deviatoricerror_groupthreshold*this->deviatoricerror_maximum;
	293	t_maxerror = this->thicknesserror_groupthreshold*this->thicknesserror_maximum;
[22241]	294	/Get the sons of the father (sibilings)/
	295	sons.clear();
	296	gmesh->Element(i)->Father()->GetHigherSubElements(sons);
	297	if(sons.size()!=4) continue;//delete just group of 4 elements. This avoids big holes in the mesh
	298	/Find the minimal distance and the error of the group/
	299	gl_groupdistance=INFINITY;if_groupdistance=INFINITY;deviatoric_grouperror=0;thickness_grouperror=0;
	300	for(int s=0;s<sons.size();s++){
	301	sid=this->index2sid[sons[s]->Index()];
	302	/Verify if this group have solutions/
	303	if(sid<0){gl_groupdistance=INFINITY;if_groupdistance=INFINITY;deviatoric_grouperror=INFINITY;thickness_grouperror=INFINITY;continue;}
	304	/Distance and error of the group/
	305	if(this->groundingline_distance>0) gl_groupdistance=std::min(gl_groupdistance,gl_distance[sid]);
	306	if(this->icefront_distance>0) if_groupdistance=std::min(if_groupdistance,if_distance[sid]);
	307	if(this->deviatoricerror_threshold>0) deviatoric_grouperror+=deviatoricerror[sid];
	308	if(this->thicknesserror_threshold>0) thickness_grouperror+=thicknesserror[sid];
	309	}
	310	/Verify the criteria/
	311	criteria=0;
[22276]	312	if(this->groundingline_distance>0 && gl_groupdistance>gl_distance_h-DBL_EPSILON) criteria++;
	313	if(this->icefront_distance>0 && if_groupdistance>if_distance_h-DBL_EPSILON) criteria++;
	314	if(this->deviatoricerror_threshold>0 && deviatoric_grouperror<d_maxerror+DBL_EPSILON) criteria++;
	315	if(this->thicknesserror_threshold>0 && thickness_grouperror<t_maxerror+DBL_EPSILON) criteria++;
[22241]	316	/Now, if the group attends the criteria, unrefine it/
	317	if(criteria==numberofcriteria){
	318	gmesh->Element(i)->Father()->ResetSubElements(); count++;
	319	for(int s=0;s<sons.size();s++){this->index2sid[sons[s]->Index()]=-1;gmesh->DeleteElement(sons[s],sons[s]->Index());}
	320	}
	321	}
	322	if(verbose) _printf_(""<<count<<"\n");
	323	/Adjust the connectivities before continue/
[23469]	324	//gmesh->BuildConnectivity(); this is not necessary
[22241]	325	/}}}/
[23066]	326
[22241]	327	/Refinement process: loop over level of refinements{{{/
	328	if(verbose) _printf_("\trefinement process (level max = "<<this->level_max<<")\n");
	329	if(verbose) _printf_("\ttotal: ");
	330	count=0;
	331	nelem=gmesh->NElements();//must keep here
	332	for(int i=0;i<nelem;i++){
	333	if(!gmesh->Element(i)) continue;
	334	if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
	335	if(gmesh->Element(i)->HasSubElement()) continue;
	336	if(gmesh->Element(i)->Level()==this->level_max) continue;
	337	/Verify if this element has solutions/
	338	sid=this->index2sid[gmesh->Element(i)->Index()];
	339	if(sid<0) continue;
[22276]	340	/Set the distance for level h/
	341	gl_distance_h = gl_distance_hmax*std::pow(this->gradation,this->level_max-(gmesh->Element(i)->Level()+1));//+1: current element level is <level_max
	342	if_distance_h = if_distance_hmax*std::pow(this->gradation,this->level_max-(gmesh->Element(i)->Level()+1));//+1: current element level is <level_max
	343	d_maxerror = this->deviatoricerror_threshold*this->deviatoricerror_maximum;
	344	t_maxerror = this->thicknesserror_threshold*this->thicknesserror_maximum;
[22241]	345	/Verify distance and error of the element, if requested/
	346	criteria=0;
[22276]	347	if(this->groundingline_distance>0 && gl_distance[sid]<gl_distance_h+DBL_EPSILON) criteria++;
	348	if(this->icefront_distance>0 && if_distance[sid]<if_distance_h+DBL_EPSILON) criteria++;
	349	if(this->deviatoricerror_threshold>0 && deviatoricerror[sid]>d_maxerror-DBL_EPSILON) criteria++;
	350	if(this->thicknesserror_threshold>0 && thicknesserror[sid]>t_maxerror-DBL_EPSILON) criteria++;
[22241]	351	/Now, if it attends any criterion, keep the element index to refine in next step/
	352	if(criteria)index.push_back(i);
	353	}
	354	/Now, refine the elements/
	355	for(int i=0;i<index.size();i++){
	356	if(!gmesh->Element(index[i])) DebugStop();
	357	if(!gmesh->Element(index[i])->HasSubElement()){gmesh->Element(index[i])->Divide(sons);count++;}
	358	}
	359	if(verbose) _printf_(""<<count<<"\n");
	360	/Adjust the connectivities before continue/
[23469]	361	//gmesh->BuildConnectivity();//this is not necessary
[22241]	362	/}}}/
[21919]	363
[22241]	364	/Now, apply smoothing and insert special elements to avoid hanging nodes{{{/
	365	this->RefineMeshWithSmoothing(verbose,gmesh);
	366	if(this->refinement_type==0) this->previousmesh=this->CreateRefPatternMesh(gmesh);//in this case, gmesh==this->fathermesh
	367	gmesh=this->previousmesh;//previous mesh is always refined to avoid hanging nodes
	368	this->RefineMeshToAvoidHangingNodes(verbose,gmesh);
	369	/}}}/
[21484]	370	}
	371	/}}}/
[25454]	372	int AmrNeopz::VerifyRefinementType(TPZGeoEl* geoel,TPZGeoMesh* gmesh){/{{{/
[22241]	373
	374	/*
	375	* 0 : no refinement
	376	* 1 : special refinement (to avoid hanging nodes)
	377	* 2 : uniform refinment
	378	* */
	379	if(!geoel) _error_("geoel is NULL!\n");
[23066]	380
[22241]	381	/Output/
	382	int type=0;
[23066]	383
[22100]	384	/Intermediaries/
	385	TPZVec<TPZGeoEl*> sons;
[23066]	386
[22241]	387	/Loop over neighboors (sides 3, 4 and 5)/
	388	for(int j=3;j<6;j++){
[23469]	389	if(!gmesh->Element(geoel->NeighbourIndex(j))->HasSubElement()) continue;
[22241]	390	sons.clear();
[23469]	391	gmesh->Element(geoel->NeighbourIndex(j))->GetHigherSubElements(sons);
	392	if(sons.size()) type++; //if neighbour was refined
	393	if(sons.size()>4) type++; //if neighbour's level is > element level+1
	394	if(type>1) break;
[22241]	395	}
[23469]	396
[22241]	397	/Verify and return/
[23469]	398	if(type>1) type=2;
[23066]	399
[22241]	400	return type;
	401	}
	402	/}}}/
[25454]	403	void AmrNeopz::RefineMeshWithSmoothing(bool &verbose,TPZGeoMesh* gmesh){/{{{/
[22100]	404
[22241]	405	/Intermediaries/
	406	int nelem =-1;
	407	int count =-1;
	408	int type =-1;
	409	int typecount =-1;
[22100]	410
[22241]	411	TPZVec<TPZGeoEl*> sons;
	412
	413	/Refinement process: loop over level of refinements/
	414	if(verbose) _printf_("\tsmoothing process (level max = "<<this->level_max<<")\n");
	415	if(verbose) _printf_("\ttotal: ");
	416
	417	count=1;
[23469]	418
[22241]	419	while(count>0){
[22100]	420	count=0;
	421	nelem=gmesh->NElements();//must keep here
[22241]	422	for(int i=0;i<nelem;i++){
[22100]	423	if(!gmesh->Element(i)) continue;
	424	if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
	425	if(gmesh->Element(i)->HasSubElement()) continue;
[22241]	426	if(gmesh->Element(i)->Level()==this->level_max) continue;
	427	/loop over neighboors (sides 3, 4 and 5). Important: neighbours has the same dimension of the element/
[23469]	428	type=this->VerifyRefinementType(gmesh->Element(i),gmesh);
[22241]	429	if(type<2){
	430	typecount=0;
	431	for(int j=3;j<6;j++){
[23469]	432	if(gmesh->Element(gmesh->Element(i)->NeighbourIndex(j))->HasSubElement()) continue;
	433	if(gmesh->Element(i)->NeighbourIndex(j)==i) typecount++;//neighbour==this element, element at the border
	434	if(this->VerifyRefinementType(gmesh->Element(gmesh->Element(i)->NeighbourIndex(j)),gmesh)==1) typecount++;
	435	if(typecount>1 && type==1) type=2;
	436	else if(typecount>2 && type==0) type=2;
	437	if(type==2) break;
[22100]	438	}
	439	}
[23469]	440
[22241]	441	/refine the element if requested/
	442	if(type==2){gmesh->Element(i)->Divide(sons); count++;}
[22100]	443	}
[22241]	444	if(verbose){
	445	if(count==0) _printf_(""<<count<<"\n");
	446	else _printf_(""<<count<<", ");
	447	}
	448	/Adjust the connectivities before continue/
[23469]	449	//gmesh->BuildConnectivity();//this is not necessary
[22100]	450	}
	451	}
	452	/}}}/
[25454]	453	void AmrNeopz::RefineMeshToAvoidHangingNodes(bool &verbose,TPZGeoMesh* gmesh){/{{{/
[23066]	454
[21919]	455	/Now, insert special elements to avoid hanging nodes/
[22100]	456	if(verbose) _printf_("\trefining to avoid hanging nodes (total: ");
[23066]	457
[22100]	458	/Intermediaries/
	459	int nelem=-1;
	460	int count= 1;
[23066]	461
[22100]	462	while(count>0){
	463	nelem=gmesh->NElements();//must keep here
	464	count=0;
	465	for(int i=0;i<nelem;i++){
	466	/Get geometric element and verify if it has already been refined. Geoel may not have been previously refined/
	467	TPZGeoEl * geoel=gmesh->Element(i);
	468	if(!geoel) continue;
	469	if(geoel->HasSubElement()) continue;
	470	if(geoel->MaterialId() != this->GetElemMaterialID()) continue;
	471	/Get the refinement pattern for this element and refine it/
	472	TPZAutoPointer<TPZRefPattern> refp=TPZRefPatternTools::PerfectMatchRefPattern(geoel);
	473	if(refp){
	474	/Non-uniform refinement/
	475	TPZVec<TPZGeoEl *> sons;
	476	geoel->SetRefPattern(refp);
	477	geoel->Divide(sons);
	478	count++;
	479	/Keep the index of the special elements/
	480	for(int j=0;j<sons.size();j++) this->specialelementsindex.push_back(sons[j]->Index());
	481	}
[21919]	482	}
[22100]	483	if(verbose){
	484	if(count==0) _printf_(""<<count<<")\n");
	485	else _printf_(""<<count<<", ");
	486	}
[23469]	487	//gmesh->BuildConnectivity();//this is not necessary
[21806]	488	}
[22100]	489	}
	490	/}}}/
[25454]	491	void AmrNeopz::DeleteSpecialElements(bool &verbose,TPZGeoMesh* gmesh){/{{{/
[22100]	492
	493	/Intermediaries/
	494	int count=0;
	495
	496	if(verbose) _printf_("\tdelete "<<this->specialelementsindex.size()<<" special elements (total: ");
	497	for(int i=0;i<this->specialelementsindex.size();i++){
	498	if(this->specialelementsindex[i]==-1) continue;
	499	if(!gmesh->Element(this->specialelementsindex[i])) continue;
	500	if(gmesh->Element(this->specialelementsindex[i])->Father()) gmesh->Element(this->specialelementsindex[i])->Father()->ResetSubElements();
	501	gmesh->DeleteElement(gmesh->Element(this->specialelementsindex[i]),this->specialelementsindex[i]);
[22241]	502	this->index2sid[this->specialelementsindex[i]]=-1;
[22100]	503	count++;
	504	}
	505	if(verbose) _printf_(""<<count<<")\n");
	506	/Cleanup/
	507	this->specialelementsindex.clear();
	508	/Adjust connectivities/
[23469]	509	//gmesh->BuildConnectivity();//this is not necessary
[21484]	510	}
	511	/}}}/
[25454]	512	void AmrNeopz::GetMesh(TPZGeoMesh* gmesh,int pdata,double pxy, int** pelements){/{{{/
[21484]	513
[21806]	514	/* IMPORTANT! pelements are in Matlab indexing
	515	NEOPZ works only in C indexing.
[21919]	516	This method cleans up and updated the this->sid2index
	517	and this->index2sid and fills in it with the new mesh.
[21806]	518	Avoid to call this method before Refinement Process.*/
[21528]	519
[21806]	520	/Intermediaries /
	521	long sid,nodeindex;
[22100]	522	int nconformelements,nconformvertices;
	523	int ntotalvertices = gmesh->NNodes();
[21806]	524	int* newelements = NULL;
[23065]	525	int* newdata = NULL;
	526	double* newmeshXY = NULL;
[21806]	527	TPZGeoEl* geoel = NULL;
[22100]	528	long* vertex_index2sid = xNew<long>(ntotalvertices);
[21919]	529	this->index2sid.clear(); this->index2sid.resize(gmesh->NElements());
[21806]	530	this->sid2index.clear();
[23066]	531
[21806]	532	/Fill in the vertex_index2sid vector with non usual index value/
	533	for(int i=0;i<gmesh->NNodes();i++) vertex_index2sid[i]=-1;
[23066]	534
[21919]	535	/Fill in the this->index2sid vector with non usual index value/
	536	for(int i=0;i<gmesh->NElements();i++) this->index2sid[i]=-1;
[23066]	537
[21806]	538	/Get elements without sons and fill in the vertex_index2sid with used vertices (indexes) /
	539	sid=0;
	540	for(int i=0;i<gmesh->NElements();i++){//over gmesh elements index
	541	geoel=gmesh->ElementVec()[i];
	542	if(!geoel) continue;
	543	if(geoel->HasSubElement()) continue;
	544	if(geoel->MaterialId() != this->GetElemMaterialID()) continue;
[22100]	545	this->sid2index.push_back(geoel->Index());//keep the element index
	546	this->index2sid[geoel->Index()]=this->sid2index.size()-1;//keep the element sid
	547	for(int j=0;j<this->GetNumberOfNodes();j++){
[21806]	548	nodeindex=geoel->NodeIndex(j);
[22100]	549	if(nodeindex<0) _error_("nodeindex is <0\n");
	550	if(vertex_index2sid[nodeindex]==-1){
[21806]	551	vertex_index2sid[nodeindex]=sid;
	552	sid++;
	553	}
	554	}
	555	}
[21484]	556
[23065]	557	/* Create new mesh structure and fill it */
[21806]	558	nconformelements = (int)this->sid2index.size();
	559	nconformvertices = (int)sid;
[22100]	560	newelements = xNew<int>(nconformelements*this->GetNumberOfNodes());
[23065]	561	newmeshXY = xNew<double>(nconformvertices*2);
	562	newdata = xNew<int>(2);
	563	newdata[0] = nconformvertices;
	564	newdata[1] = nconformelements;
[21806]	565
[22100]	566	for(int i=0;i<ntotalvertices;i++){//over the TPZNode index (fill in the ISSM vertices coords)
	567	sid=vertex_index2sid[i];
	568	if(sid==-1) continue;//skip this index (node no used)
	569	TPZVec<REAL> coords(3,0.);
	570	gmesh->NodeVec()[i].GetCoordinates(coords);
[23065]	571	newmeshXY[2*sid] = coords[0]; // X
	572	newmeshXY[2*sid+1] = coords[1]; // Y
[21806]	573	}
[23066]	574
[21806]	575	for(int i=0;i<this->sid2index.size();i++){//over the sid (fill the ISSM elements)
[22100]	576	for(int j=0;j<this->GetNumberOfNodes();j++) {
[21806]	577	geoel = gmesh->ElementVec()[this->sid2index[i]];
	578	sid = vertex_index2sid[geoel->NodeIndex(j)];
[22100]	579	newelements[i*this->GetNumberOfNodes()+j]=(int)sid+1;//C to Matlab indexing
[21806]	580	}
[21919]	581	/*Verify the Jacobian determinant. If detJ<0, swap the 2 first postions:
	582	a -> b
	583	b -> a */
	584	double detJ,xa,xb,xc,ya,yb,yc;
	585	int a,b,c;
	586
[22100]	587	a=newelements[i*this->GetNumberOfNodes()+0]-1;
	588	b=newelements[i*this->GetNumberOfNodes()+1]-1;
	589	c=newelements[i*this->GetNumberOfNodes()+2]-1;
[21919]	590
[23065]	591	xa=newmeshXY[2a]; ya=newmeshXY[2a+1];
	592	xb=newmeshXY[2b]; yb=newmeshXY[2b+1];
	593	xc=newmeshXY[2c]; yc=newmeshXY[2c+1];
[21919]	594
	595	detJ=(xb-xa)(yc-ya)-(xc-xa)(yb-ya);
[23066]	596
[21919]	597	/verify and swap, if necessary/
	598	if(detJ<0) {
[22100]	599	newelements[i*this->GetNumberOfNodes()+0]=b+1;//a->b
	600	newelements[i*this->GetNumberOfNodes()+1]=a+1;//b->a
[21919]	601	}
[21806]	602	}
[21919]	603
[21806]	604	/Setting outputs/
[23065]	605	*pdata = newdata;
	606	*pxy = newmeshXY;
[21806]	607	*pelements = newelements;
[23066]	608
[21806]	609	/Cleanup/
	610	xDelete<long>(vertex_index2sid);
[21484]	611	}
	612	/}}}/
[25454]	613	void AmrNeopz::Initialize(){/{{{/
[21484]	614
[21806]	615	/* IMPORTANT! elements come in Matlab indexing
	616	NEOPZ works only in C indexing*/
[23066]	617
[23501]	618	if(this->numberofvertices<=0) _error_("Impossible to create initial mesh: nvertices is <= 0!\n");
	619	if(this->numberofelements<=0) _error_("Impossible to create initial mesh: nelements is <= 0!\n");
[22241]	620	if(this->refinement_type!=0 && this->refinement_type!=1) _error_("Impossible to create initial mesh: refinement type is not defined!\n");
[21484]	621
	622	/Verify and creating initial mesh/
[23501]	623	if(!this->x \|\| !this->y \|\| !this->elementslist) _error_("Mesh data structure is NULL!\n");
	624	if(this->fathermesh \|\| this->previousmesh) _error_("Initial mesh already exists!\n");
[23066]	625
[21503]	626	this->fathermesh = new TPZGeoMesh();
[23501]	627	this->fathermesh->NodeVec().Resize(this->numberofvertices);
[21484]	628
[21806]	629	/Set the vertices (geometric nodes in NeoPZ context)/
[23501]	630	for(int i=0;i<this->numberofvertices;i++){
[21528]	631	/x,y,z coords/
	632	TPZManVector<REAL,3> coord(3,0.);
[23501]	633	coord[0]= this->x[i];
	634	coord[1]= this->y[i];
[21806]	635	coord[2]= 0.;
[21528]	636	/Insert in the mesh/
	637	this->fathermesh->NodeVec()[i].SetCoord(coord);
	638	this->fathermesh->NodeVec()[i].SetNodeId(i);
[21484]	639	}
[23066]	640
[21484]	641	/Generate the elements/
[24491]	642	int64_t index;
[21528]	643	const int mat = this->GetElemMaterialID();
[24491]	644	TPZManVector<int64_t> elem(this->GetNumberOfNodes(),0);
[23501]	645	this->index2sid.clear(); this->index2sid.resize(this->numberofelements);
[21806]	646	this->sid2index.clear();
[21484]	647
[23501]	648	for(int i=0;i<this->numberofelements;i++){
	649	for(int j=0;j<this->GetNumberOfNodes();j++) elem[j]=this->elementslist[i*this->GetNumberOfNodes()+j]-1;//Convert Matlab to C indexing
[22100]	650	switch(this->GetNumberOfNodes()){
[22241]	651	case 3: this->fathermesh->CreateGeoElement(ETriangle,elem,mat,index,this->refinement_type); break;
[21806]	652	default: _error_("mesh not supported yet");
[21484]	653	}
[21528]	654	/Define the element ID/
[21806]	655	this->fathermesh->ElementVec()[index]->SetId(i);
[22241]	656	/Initialize sid2index and index2sid/
[21806]	657	this->sid2index.push_back((int)index);
[22241]	658	this->index2sid[(int)index]=this->sid2index.size()-1;//keep the element sid
[21484]	659	}
[21528]	660	/Build element and node connectivities/
	661	this->fathermesh->BuildConnectivity();
[22100]	662	/Set previous mesh/
[22241]	663	if(this->refinement_type==1) this->previousmesh=new TPZGeoMesh(*this->fathermesh);
	664	else this->previousmesh=this->CreateRefPatternMesh(this->fathermesh);
[21484]	665	}
	666	/}}}/
[25454]	667	TPZGeoMesh* AmrNeopz::CreateRefPatternMesh(TPZGeoMesh* gmesh){/{{{/
[23066]	668
[21672]	669	TPZGeoMesh *newgmesh = new TPZGeoMesh();
	670	newgmesh->CleanUp();
[23066]	671
[21672]	672	int nnodes = gmesh->NNodes();
	673	int nelem = gmesh->NElements();
	674	int mat = this->GetElemMaterialID();;
	675	int reftype = 1;
[24491]	676	int64_t index;
[22241]	677
[21672]	678	//nodes
	679	newgmesh->NodeVec().Resize(nnodes);
	680	for(int i=0;i<nnodes;i++) newgmesh->NodeVec()[i] = gmesh->NodeVec()[i];
[23066]	681
[21672]	682	//elements
	683	for(int i=0;i<nelem;i++){
	684	TPZGeoEl * geoel = gmesh->Element(i);
[23066]	685
[22241]	686	if(!geoel){
	687	index=newgmesh->ElementVec().AllocateNewElement();
	688	newgmesh->ElementVec()[index] = NULL;
	689	continue;
	690	}
[23066]	691
[24491]	692	TPZManVector<int64_t> elem(3,0);
[21672]	693	for(int j=0;j<3;j++) elem[j] = geoel->NodeIndex(j);
[23066]	694
[21672]	695	newgmesh->CreateGeoElement(ETriangle,elem,mat,index,reftype);
[22241]	696	newgmesh->ElementVec()[index]->SetId(geoel->Id());
[23066]	697
[21672]	698	TPZGeoElRefPattern<TPZGeoTriangle>* newgeoel = dynamic_cast<TPZGeoElRefPattern<TPZGeoTriangle>*>(newgmesh->ElementVec()[index]);
[23066]	699
[21672]	700	//old neighbourhood
	701	const int nsides = TPZGeoTriangle::NSides;
	702	TPZVec< std::vector<TPZGeoElSide> > neighbourhood(nsides);
	703	TPZVec<long> NodesSequence(0);
	704	for(int s = 0; s < nsides; s++){
	705	neighbourhood[s].resize(0);
	706	TPZGeoElSide mySide(geoel,s);
	707	TPZGeoElSide neighS = mySide.Neighbour();
	708	if(mySide.Dimension() == 0){
	709	long oldSz = NodesSequence.NElements();
	710	NodesSequence.resize(oldSz+1);
	711	NodesSequence[oldSz] = geoel->NodeIndex(s);
	712	}
	713	while(mySide != neighS){
	714	neighbourhood[s].push_back(neighS);
	715	neighS = neighS.Neighbour();
	716	}
	717	}
[23066]	718
[21672]	719	//inserting in new element
	720	for(int s = 0; s < nsides; s++){
	721	TPZGeoEl * tempEl = newgeoel;
	722	TPZGeoElSide tempSide(newgeoel,s);
	723	int byside = s;
	724	for(unsigned long n = 0; n < neighbourhood[s].size(); n++){
	725	TPZGeoElSide neighS = neighbourhood[s][n];
	726	tempEl->SetNeighbour(byside, neighS);
	727	tempEl = neighS.Element();
	728	byside = neighS.Side();
	729	}
	730	tempEl->SetNeighbour(byside, tempSide);
	731	}
[23066]	732
[21672]	733	long fatherindex = geoel->FatherIndex();
	734	if(fatherindex>-1) newgeoel->SetFather(fatherindex);
[23066]	735
[21672]	736	if(!geoel->HasSubElement()) continue;
[23066]	737
[21672]	738	int nsons = geoel->NSubElements();
	739
	740	TPZAutoPointer<TPZRefPattern> ref = gRefDBase.GetUniformRefPattern(ETriangle);
	741	newgeoel->SetRefPattern(ref);
[23066]	742
[21672]	743	for(int j=0;j<nsons;j++){
	744	TPZGeoEl* son = geoel->SubElement(j);
	745	if(!son){
	746	DebugStop();
	747	}
	748	newgeoel->SetSubElement(j,son);
	749	}
	750	}
[22241]	751
	752	/Now, build connectivities/
[21672]	753	newgmesh->BuildConnectivity();
[23066]	754
[21672]	755	return newgmesh;
[21484]	756	}
	757	/}}}/
[25454]	758	void AmrNeopz::CheckMesh(int pdata,double pxy,int** pelements){/{{{/
[21528]	759
[21806]	760	/Basic verification/
[23065]	761	if(!pdata) _error_("Impossible to continue: pdata is NULL!\n");
[24492]	762	if(**pdata<=0) _error_("Impossible to continue: nvertices <=0!\n");
	763	if((pdata+1)<=0) _error_("Impossible to continue: nelements <=0!\n");
[23065]	764	if(!pxy) _error_("Impossible to continue: pxy is NULL!\n");
[21806]	765	if(!pelements) _error_("Impossible to continue: pelements is NULL!\n");
[21484]	766	}
	767	/}}}/
[25454]	768	void AmrNeopz::PrintGMeshVTK(TPZGeoMesh* gmesh,std::ofstream &file,bool matColor){/{{{/
[23066]	769
[22100]	770	file.clear();
	771	long nelements = gmesh->NElements();
	772	TPZGeoEl *gel;
	773	std::stringstream node, connectivity, type, material;
	774
	775	//Header
	776	file << "# vtk DataFile Version 3.0" << std::endl;
	777	file << "TPZGeoMesh VTK Visualization" << std::endl;
	778	file << "ASCII" << std::endl << std::endl;
	779	file << "DATASET UNSTRUCTURED_GRID" << std::endl;
	780	file << "POINTS ";
	781
	782	long actualNode = -1, size = 0, nVALIDelements = 0;
	783	for(long el = 0; el < nelements; el++){
	784	gel = gmesh->ElementVec()[el];
	785	if(!gel )//\|\| (gel->Type() == EOned && !gel->IsLinearMapping()))//Exclude Arc3D and Ellipse3D
	786	{
	787	continue;
	788	}
	789	if(gel->HasSubElement())
	790	{
	791	continue;
	792	}
	793	MElementType elt = gel->Type();
	794	int elNnodes = MElementType_NNodes(elt);
[23066]	795
[22100]	796	size += (1+elNnodes);
	797	connectivity << elNnodes;
[23066]	798
[22100]	799	for(int t = 0; t < elNnodes; t++)
	800	{
	801	for(int c = 0; c < 3; c++)
	802	{
	803	double coord = gmesh->NodeVec()[gel->NodeIndex(t)].Coord(c);
	804	node << coord << " ";
	805	}
	806	node << std::endl;
[23066]	807
[22100]	808	actualNode++;
	809	connectivity << " " << actualNode;
	810	}
	811	connectivity << std::endl;
[23066]	812
[22100]	813	int elType = this->GetVTK_ElType(gel);
	814	type << elType << std::endl;
[23066]	815
[22100]	816	if(matColor == true)
	817	{
	818	material << gel->MaterialId() << std::endl;
	819	}
	820	else
	821	{
	822	material << gel->Index() << std::endl;
	823	}
[23066]	824
[22100]	825	nVALIDelements++;
	826	}
	827	node << std::endl;
	828	actualNode++;
	829	file << actualNode << " float" << std::endl << node.str();
	830
	831	file << "CELLS " << nVALIDelements << " ";
	832
	833	file << size << std::endl;
	834	file << connectivity.str() << std::endl;
	835
	836	file << "CELL_TYPES " << nVALIDelements << std::endl;
	837	file << type.str() << std::endl;
	838
	839	file << "CELL_DATA" << " " << nVALIDelements << std::endl;
	840	file << "FIELD FieldData 1" << std::endl;
	841	if(matColor == true)
	842	{
	843	file << "material 1 " << nVALIDelements << " int" << std::endl;
	844	}
	845	else
	846	{
	847	file << "ElementIndex 1 " << nVALIDelements << " int" << std::endl;
	848	}
	849	file << material.str();
	850	file.close();
	851	}
	852	/}}}/
[25454]	853	int AmrNeopz::GetVTK_ElType(TPZGeoEl * gel){/{{{/
[23066]	854
[22100]	855	MElementType pzElType = gel->Type();
[23066]	856
[22100]	857	int elType = -1;
	858	switch (pzElType)
	859	{
	860	case(EPoint):
	861	{
	862	elType = 1;
	863	break;
	864	}
	865	case(EOned):
	866	{
	867	elType = 3;
	868	break;
	869	}
	870	case (ETriangle):
	871	{
	872	elType = 5;
	873	break;
	874	}
	875	case (EQuadrilateral):
	876	{
	877	elType = 9;
	878	break;
	879	}
	880	case (ETetraedro):
	881	{
	882	elType = 10;
	883	break;
	884	}
	885	case (EPiramide):
	886	{
	887	elType = 14;
	888	break;
	889	}
	890	case (EPrisma):
	891	{
	892	elType = 13;
	893	break;
	894	}
	895	case (ECube):
	896	{
	897	elType = 12;
	898	break;
	899	}
	900	default:
	901	{
	902	std::cout << "Element type not found on " << __PRETTY_FUNCTION__ << std::endl;
	903	DebugStop();
	904	break;
	905	}
	906	}
	907	if(elType == -1)
	908	{
	909	std::cout << "Element type not found on " << __PRETTY_FUNCTION__ << std::endl;
	910	std::cout << "MIGHT BE CURVED ELEMENT (quadratic or quarter point)" << std::endl;
	911	DebugStop();
	912	}
[23066]	913
[22100]	914	return elType;
	915	}
	916	/}}}/
[25454]	917	void AmrNeopz::ReadMesh(){/{{{/
[23471]	918
	919	std::string fathermeshfile = "/home/santos/issm_fathermesh.txt";
	920	std::string previousmeshfile = "/home/santos/issm_previousmesh.txt";
	921	std::string amrfile = "/home/santos/issm_amr.txt";
	922	std::ifstream amrifstream(amrfile.c_str());
	923	int size,value;
[23501]	924
[23471]	925	TPZPersistenceManager::OpenRead(fathermeshfile);
	926	this->fathermesh = dynamic_cast<TPZGeoMesh *>(TPZPersistenceManager::ReadFromFile());
	927	TPZPersistenceManager::CloseRead();
	928
	929	TPZPersistenceManager::OpenRead(previousmeshfile);
	930	this->previousmesh = dynamic_cast<TPZGeoMesh *>(TPZPersistenceManager::ReadFromFile());
	931	TPZPersistenceManager::CloseRead();
	932
	933	if(!amrifstream.is_open()) _error_("amr ifstream is not open!");
	934	amrifstream.seekg(0);
	935
	936	this->sid2index.clear();
	937	this->index2sid.clear();
	938	this->specialelementsindex.clear();
	939
	940	amrifstream>>size;
	941	for(int i=0;i<size;i++){
	942	amrifstream>>value;
	943	this->sid2index.push_back(value);
	944	}
	945
	946	amrifstream>>size;
	947	for(int i=0;i<size;i++){
	948	amrifstream>>value;
	949	this->index2sid.push_back(value);
	950	}
	951
	952	amrifstream>>size;
	953	for(int i=0;i<size;i++){
	954	amrifstream>>value;
	955	this->specialelementsindex.push_back(value);
	956	}
	957	}
	958	/}}}/
[25454]	959	void AmrNeopz::WriteMesh(){/{{{/
[23471]	960
	961	std::string fathermeshfile = "/home/santos/issm_fathermesh.txt";
	962	std::string previousmeshfile = "/home/santos/issm_previousmesh.txt";
	963	std::string amrfile = "/home/santos/issm_amr.txt";
	964	std::ofstream amrofstream(amrfile.c_str());
	965
	966	if(this->fathermesh){
	967	TPZPersistenceManager::OpenWrite(fathermeshfile);
	968	TPZPersistenceManager::WriteToFile(this->fathermesh);
	969	TPZPersistenceManager::CloseWrite();
	970	}
	971
	972	if(this->previousmesh){
	973	TPZPersistenceManager::OpenWrite(previousmeshfile);
	974	TPZPersistenceManager::WriteToFile(this->previousmesh);
	975	TPZPersistenceManager::CloseWrite();
	976	}
	977
	978	if(this->sid2index.size()>0){
	979	amrofstream << this->sid2index.size() << std::endl;
	980	for(int i=0;i<this->sid2index.size();i++) {
	981	amrofstream << this->sid2index[i] << std::endl;
	982	}
	983	}
	984
	985	if(this->index2sid.size()>0){
	986	amrofstream << this->index2sid.size() << std::endl;
	987	for(int i=0;i<this->index2sid.size();i++) {
	988	amrofstream << this->index2sid[i] << std::endl;
	989	}
	990	}
	991
	992	if(this->specialelementsindex.size()){
	993	amrofstream << this->specialelementsindex.size() << std::endl;
	994	for(int i=0;i<this->specialelementsindex.size();i++) {
	995	amrofstream << this->specialelementsindex[i] << std::endl;
	996	}
	997	}
	998	amrofstream.flush();
	999	amrofstream.close();
	1000	}
	1001	/}}}/

Note: See TracBrowser for help on using the repository browser.

Download in other formats: