Context Navigation

source: issm/trunk/src/c/classes/AdaptiveMeshRefinement.cpp@ 23189

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

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

Download in other formats: