Index: issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp	(revision 21806)
+++ issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp	(revision 21919)
@@ -10,8 +10,5 @@
 
 #include "./AdaptiveMeshRefinement.h"
-#include "TPZVTKGeoMesh.h"
-#include "../shared/shared.h"
-#include "pzgeotriangle.h"
-#include "pzreftriangle.h"
+
 using namespace pzgeom;
 
@@ -40,23 +37,15 @@
 	this->sid2index.resize(cp.sid2index.size());
 	for(int i=0;i<cp.sid2index.size();i++) this->sid2index[i]=cp.sid2index[i];
-	
+	this->index2sid.clear();
+	this->index2sid.resize(cp.index2sid.size());
+	for(int i=0;i<cp.index2sid.size();i++) this->index2sid[i]=cp.index2sid[i];
+	this->specialelementsindex.clear();
+	this->specialelementsindex.resize(cp.specialelementsindex.size());
+	for(int i=0;i<cp.specialelementsindex.size();i++) this->specialelementsindex[i]=cp.specialelementsindex[i];
+
 	return *this;
-
 }
 /*}}}*/
 AdaptiveMeshRefinement::~AdaptiveMeshRefinement(){/*{{{*/
-	
-	bool ismismip = false;
-	if(ismismip){//itapopo
-		TPZFileStream fstr;
-		std::stringstream ss;
-	    
-		ss << this->levelmax;
-		std::string AMRfile	= "/home/santos/L" + ss.str() + "_amr.txt"; 
-	
-		fstr.OpenWrite(AMRfile.c_str());
-		int withclassid = 1;
-		this->Write(fstr,withclassid);
-	}
 	this->CleanUp();
 	gRefDBase.clear();
@@ -73,4 +62,6 @@
 	this->regionlevelmax		= -1;
 	this->sid2index.clear();
+	this->index2sid.clear();
+	this->specialelementsindex.clear();
 }
 /*}}}*/
@@ -84,79 +75,9 @@
 	this->regionlevel1		= -1;
 	this->regionlevelmax		= -1;
+	this->step					= 0;
+	this->smooth_frequency	= 1;
 	this->sid2index.clear();
-}
-/*}}}*/
-int AdaptiveMeshRefinement::ClassId() const{/*{{{*/
-    return 13829430; //Antartic area with ice shelves (km^2)
-}
-/*}}}*/
-void AdaptiveMeshRefinement::Read(TPZStream &buf,void *context){/*{{{*/
-
-	try
-	{
-		/* Read the id context*/
-		TPZSaveable::Read(buf,context);
-		/* Read class id*/
-		int classid;
-		buf.Read(&classid,1); 
-		/* Verify the class id*/
-      if(classid!=this->ClassId()) _error_("AdaptiveMeshRefinement::Read: Error in restoring AdaptiveMeshRefinement!\n"); 
-		/* Read simple attributes */
-		buf.Read(&this->levelmax,1);
-		buf.Read(&this->elementswidth,1);
-		buf.Read(&this->regionlevel1,1);
-		buf.Read(&this->regionlevelmax,1);
-		/* Read vector attributes*/
-		int size;
-		buf.Read(&size,1);
-		int* psid2index=xNew<int>(size);
-		buf.Read(psid2index,size);
-		this->sid2index.clear();
-		this->sid2index.assign(psid2index,psid2index+size);
-		/* Read geometric mesh (father)*/
-		TPZSaveable *sv1 = TPZSaveable::Restore(buf,0);
-		this->fathermesh = dynamic_cast<TPZGeoMesh*>(sv1);
-		/* Read geometric mesh (current)*/
-		TPZSaveable *sv2 = TPZSaveable::Restore(buf,0);
-		this->currentmesh = dynamic_cast<TPZGeoMesh*>(sv2);
-		/* Cleanup*/
-		xDelete<int>(psid2index);
-	}
-	catch(const std::exception& e)
-	{
-		_error_("AdaptiveMeshRefinement::Read: Exception catched!\n");
-	}
-
-}
-/*}}}*/
-template class TPZRestoreClass<AdaptiveMeshRefinement,13829430>;/*{{{*/
-/*}}}*/
-void AdaptiveMeshRefinement::Write(TPZStream &buf,int withclassid){/*{{{*/
-    
-	try
-	{
-		/* Write context (this class) class ID*/
-		TPZSaveable::Write(buf,withclassid);
-		/* Write this class id*/
-		int classid = this->ClassId();
-		buf.Write(&classid,1);
-		/* Write simple attributes */
-		buf.Write(&this->levelmax,1);
-		buf.Write(&this->elementswidth,1);
-		buf.Write(&this->regionlevel1,1);
-		buf.Write(&this->regionlevelmax,1);
-		/* Write vector attributes*/
-		int size=this->sid2index.size();
-		int* psid2index=&this->sid2index[0];
-		buf.Write(&size,1);//vector size
-		buf.Write(psid2index,this->sid2index.size());
-		/* Write the geometric mesh*/
-		this->fathermesh->Write(buf,this->ClassId());
-		this->currentmesh->Write(buf,this->ClassId());
-    }
-    catch(const std::exception& e)
-    {
-		_error_("AdaptiveMeshRefinement::Write: Exception catched!\n");
-    }
+	this->index2sid.clear();
+	this->specialelementsindex.clear();
 }
 /*}}}*/
@@ -176,27 +97,12 @@
 	if(numberofelements!=this->sid2index.size()) _error_("Impossible to execute refinement: sid2index.size is not equal to numberofelements!\n");
 
-	/*Execute the refinement.*/
-	this->RefinementProcess(amr_verbose,partiallyfloatedelements,masklevelset,deviatorictensorerror,thicknesserror);
-    
-	/*Get new geometric mesh in ISSM data structure*/
-	this->GetMesh(newnumberofvertices,newnumberofelements,x,y,elementslist);
-	
-	/*Verify the new geometry*/
-	this->CheckMesh(newnumberofvertices,newnumberofelements,this->elementswidth,x,y,elementslist);
-
-}
-/*}}}*/
-void AdaptiveMeshRefinement::RefinementProcess(bool &amr_verbose,double* partiallyfloatedelements,double* masklevelset,
-																double* deviatorictensorerror,double* thicknesserror){/*{{{*/
-   
-	if(amr_verbose) _printf_("\n\trefinement process started (level max = " << this->levelmax << ")\n");
-	
-	/*Intermediaries*/
-	TPZGeoMesh* nohangingnodesmesh=NULL;
+	/*Combine the fields*/
 	double mean_mask		= 0;
 	double mean_tauerror = 0;
 	double mean_Herror	= 0;
-	double group_error	= 0;
-   
+	int* fmask				= xNew<int>(numberofelements);
+	int* ftauerror			= xNew<int>(numberofelements);
+	int* fHerror			= xNew<int>(numberofelements);
+	int* phi					= xNew<int>(numberofelements);
 	/*Calculate mean values*/
 	for(int i=0;i<this->sid2index.size();i++){
@@ -208,32 +114,183 @@
 	mean_tauerror	/= this->sid2index.size();
 	mean_Herror		/= this->sid2index.size();
-
+	/*Filter to thickness*/
+	for(int i=0;i<this->sid2index.size();i++){
+		fmask[i]=0;
+		if(thicknesserror[i]>mean_Herror) fmask[i]=1;
+	}
+	/*Filter to tau*/
+	for(int i=0;i<this->sid2index.size();i++){
+		ftauerror[i]=0;
+		if(deviatorictensorerror[i]>mean_tauerror) ftauerror[i]=1;
+	}
+	/*Sum*/
+	for(int i=0;i<this->sid2index.size();i++){
+		phi[i]=ftauerror[i]+fmask[i];
+	}
+
+	/*Execute the refinement.*/
+	this->RefinementProcess(amr_verbose,partiallyfloatedelements,masklevelset,deviatorictensorerror,thicknesserror);
+   
+	/*Get new geometric mesh in ISSM data structure*/
+	this->GetMesh(newnumberofvertices,newnumberofelements,x,y,elementslist);
+
+	std::ofstream file("/home/santos/mesh.vtk");
+	TPZVTKGeoMesh::PrintGMeshVTK(this->currentmesh,file);
+
+	/*Verify the new geometry*/
+	this->CheckMesh(newnumberofvertices,newnumberofelements,this->elementswidth,x,y,elementslist);
+
+}
+/*}}}*/
+void AdaptiveMeshRefinement::RefinementProcess(bool &amr_verbose,double* partiallyfloatedelements,double* masklevelset,double* deviatorictensorerror,double* thicknesserror){/*{{{*/
+   
+	if(amr_verbose) _printf_("\n\trefinement process started (level max = " << this->levelmax << ")\n");
+	
+	/*Intermediaries*/
+	TPZGeoMesh* nohangingnodesmesh=NULL;
+
+	//itapopo
+	this->step++;
+
+	double mean_mask		= 0;
+	double mean_tauerror = 0;
+	double mean_Herror	= 0;
+	double group_Herror	= 0;
+	int index,sid;
+	std::vector<int> specialfatherindex; specialfatherindex.clear();
+
+	/*Calculate mean values*/
+	for(int i=0;i<this->sid2index.size();i++){
+		mean_mask		+= masklevelset[i]; 
+		mean_tauerror	+= deviatorictensorerror[i]; 
+		mean_Herror		+= thicknesserror[i];
+	}
+	mean_mask		/= this->sid2index.size();
+	mean_tauerror	/= this->sid2index.size();
+	mean_Herror		/= this->sid2index.size();
+
+	if(amr_verbose) _printf_("\t\tdeal with special elements...\n");
+	/*Deal with special elements*/
+	for(int i=0;i<this->specialelementsindex.size();i++){
+		if(this->specialelementsindex[i]==-1) continue;
+		/*Get special element and verify*/
+		TPZGeoEl* geoel=this->currentmesh->Element(this->specialelementsindex[i]);
+		if(!geoel)_error_("special element (sid) "<<i<<" is null!\n");
+		if(geoel->HasSubElement())_error_("special element (sid) "<<i<<" has "<<geoel->NSubElements()<<" subelements!\n");
+		if(geoel->MaterialId()!=this->GetElemMaterialID()) _error_("geoel->MaterialId is not GetElemMaterialID!\n");
+		if(!geoel->Father())_error_("father of special element (sid) "<<i<<" is null!\n");
+		
+		/*Get element's siblings and verify*/
+		TPZGeoEl* father=geoel->Father();
+		TPZVec<TPZGeoEl *> siblings;
+		father->GetHigherSubElements(siblings);
+		std::vector<int> sidvec; sidvec.resize(siblings.size());
+		for (int j=0;j<siblings.size();j++){
+			if(!siblings[j]) _error_("special element (sid) "<<i<<" has a null siblings null!\n"); 
+			sidvec[j]=this->index2sid[siblings[j]->Index()];
+		}
+		
+		/*Now, reset the data strucure and verify if the siblings should be deleted*/	
+		if(siblings.size()<4){
+			/*Reset subelements in the father*/
+			father->ResetSubElements();
+		}else{
+			if(siblings.size()!=4) _error_("element (index) "<<father->Index()<<" has "<<father->NSubElements()<<" subelements!\n");
+		}
+		for (int j=0;j<siblings.size();j++){
+			for(int k=0;k<this->specialelementsindex.size();k++){
+				if(this->specialelementsindex[k]==siblings[j]->Index()){ 
+					index									= siblings[j]->Index();
+					if(index<0) _error_("index is null!\n");
+					sid									= this->index2sid[index];
+					if(sid<0) _error_("sid is null!\n");
+					this->specialelementsindex[k]	= -1;
+					this->index2sid[index]			= -1;
+					this->sid2index[sid]				= -1;
+				}
+			}
+			if(siblings.size()<4){
+				/*Ok, the special element can be deleted*/
+				siblings[j]->ResetSubElements();
+				this->currentmesh->DeleteElement(siblings[j],siblings[j]->Index());
+			}
+		}
+		
+		/*Now, verify if the father should be refined with uniform pattern (smoother)*/
+		this->smooth_frequency=10000000;
+		if(siblings.size()==3 || this->step%this->smooth_frequency==0){//it keeps the mesh with uniform elements
+			/*Father has uniform subelements now*/
+			TPZVec<TPZGeoEl *> sons;
+			father->Divide(sons);
+			this->smooth_frequency=this->step;
+		}else{
+			specialfatherindex.push_back(father->Index());
+		}
+		if(this->specialelementsindex[i]!=-1) _error_("special element "<<i<<" was not deleted!\n");	
+	}
+	this->currentmesh->BuildConnectivity();
+	
 	if(amr_verbose) _printf_("\t\tuniform refinement...\n");
+	/*Deal with uniform elemnts*/
 	for(int i=0;i<this->sid2index.size();i++){
+		if(this->sid2index[i]==-1) continue;
+		/*Get element and verify*/
 		TPZGeoEl* geoel=this->currentmesh->Element(this->sid2index[i]);
-		if(geoel->HasSubElement()) _error_("Impossible to refine: geoel (index) "<<this->sid2index[i]<<" has subelements!\n");
-		if(geoel->MaterialId()!=this->GetElemMaterialID()) _error_("Impossible to refine: geoel->MaterialId is not GetElemMaterialID!\n");
-		/*Refine*/
-		if(thicknesserror[i]>mean_Herror){
+		if(geoel->HasSubElement()) _error_("element (sid) "<<i<<" has "<<geoel->NSubElements()<<" subelements!\n");
+		if(geoel->MaterialId()!=this->GetElemMaterialID()) _error_("geoel->MaterialId is not GetElemMaterialID!\n");
+
+		/*Refine process*/
+		if(thicknesserror[i]>mean_Herror)
+		{	
+			int count=0;
+			TPZGeoEl* father=geoel->Father();
+			if(father){
+				for(int j=3;j<6;j++){
+					index=father->Neighbour(j).Element()->Index();
+					for(int k=0;k<specialfatherindex.size();k++) if(specialfatherindex[k]==index) count++;
+				}
+			}
 			TPZVec<TPZGeoEl *> sons;
-			if(geoel->Level()<this->levelmax) geoel->Divide(sons);
+			if(geoel->Level()<this->levelmax && count==0) geoel->Divide(sons);
 		} 
-		else if(geoel->Level()>0){ /*try to unrefine*/
-			TPZVec<TPZGeoEl *> sons;
-			geoel->Father()->GetHigherSubElements(sons);
-			group_error=0;
-			for(int j=0;j<sons.size();j++){
-				sons[j]->Index();
+		else if(geoel->Level()>0)
+		{/*Unrefine process*/
+			
+			/*Get siblings and verify*/
+			TPZVec<TPZGeoEl *> siblings;
+			geoel->Father()->GetHigherSubElements(siblings);
+			//if(siblings.size()<4) _error_("Impossible to refine: geoel (index) "<<this->sid2index[i]<<" has less than 3 siblings!\n");	
+			if(siblings.size()>4) continue;//Father has more then 4 sons, this group should not be unrefined.
+			
+			/*Compute the error of the group*/
+			group_Herror=0;
+			for(int j=0;j<siblings.size();j++){
+				index		= siblings[j]->Index();
+				sid		= this->index2sid[index];
+				if(sid==-1) continue;
+				group_Herror+=thicknesserror[sid];
 			}
-		}
-	}
+			/*Verify if this group should be unrefined*/
+			if(group_Herror>0 && group_Herror<0*mean_Herror){ //itapopo
+				/*Reset subelements in the father*/
+				this->currentmesh->Element(geoel->Father()->Index())->ResetSubElements();
+				/*Delete the elements and set their indexes in the index2sid and sid2index*/
+				for (int j=0;j<siblings.size();j++){
+					index	= siblings[j]->Index();
+					sid	= this->index2sid[index];
+					this->index2sid[index]=-1;
+					if(sid!=-1) this->sid2index[sid]=-1;
+					this->currentmesh->DeleteElement(siblings[j],siblings[j]->Index());
+				}//for j
+			}//if
+		}/*Unrefine process*/
+	}//for i
 	this->currentmesh->BuildConnectivity();
 	
 	if(amr_verbose) _printf_("\t\trefine to avoid hanging nodes...\n");
 	this->RefineMeshToAvoidHangingNodes(this->currentmesh);
-	this->currentmesh->BuildConnectivity();
 	
 		//nohangingnodesmesh = this->CreateRefPatternMesh(newmesh); itapopo tentar otimizar
-	
+
 	if(amr_verbose) _printf_("\trefinement process done!\n");
 }
@@ -257,5 +314,6 @@
 void AdaptiveMeshRefinement::RefineMeshToAvoidHangingNodes(TPZGeoMesh *gmesh){/*{{{*/
    
-	/*Refine elements to avoid hanging nodes: non-uniform refinement*/
+	/*Now, insert special elements to avoid hanging nodes*/
+	this->specialelementsindex.clear();
 	const int NElem = gmesh->NElements();
 	for(int i=0;i<NElem;i++){
@@ -268,11 +326,13 @@
 		TPZAutoPointer<TPZRefPattern> refp=TPZRefPatternTools::PerfectMatchRefPattern(geoel);
 		if(refp){
+			/*Non-uniform refinement*/
 			TPZVec<TPZGeoEl *> Sons;
 			geoel->SetRefPattern(refp);
 			geoel->Divide(Sons);
-      }
-	}
-   gmesh->BuildConnectivity();
-    
+			/*Keep the index of the special elements*/
+			for(int j=0;j<Sons.size();j++) this->specialelementsindex.push_back(Sons[j]->Index());
+		}
+	}
+	gmesh->BuildConnectivity();
 }
 /*}}}*/
@@ -281,6 +341,6 @@
 	/* IMPORTANT! pelements are in Matlab indexing
 	   NEOPZ works only in C indexing.
-		This method cleans up and updated the this->sid2index 
-		and fills in it with the new mesh.
+		This method cleans up and updated the this->sid2index
+		and this->index2sid and fills in it with the new mesh.
 		Avoid to call this method before Refinement Process.*/
 
@@ -295,4 +355,5 @@
 	TPZGeoEl* geoel			= NULL;
 	long* vertex_index2sid 	= xNew<long>(gmesh->NNodes());
+	this->index2sid.clear(); this->index2sid.resize(gmesh->NElements());
 	this->sid2index.clear();
 	
@@ -307,4 +368,7 @@
 	/*Fill in the vertex_index2sid vector with non usual index value*/
 	for(int i=0;i<gmesh->NNodes();i++) vertex_index2sid[i]=-1;
+	
+	/*Fill in the this->index2sid vector with non usual index value*/
+	for(int i=0;i<gmesh->NElements();i++) this->index2sid[i]=-1;
 	
 	/*Get elements without sons and fill in the vertex_index2sid with used vertices (indexes) */
@@ -316,4 +380,5 @@
 		if(geoel->MaterialId() != this->GetElemMaterialID()) continue;
 		this->sid2index.push_back(i);//keep the element index
+		this->index2sid[i]=this->sid2index.size()-1;//keep the element sid
 		for(int j=0;j<this->elementswidth;j++){
       	nodeindex=geoel->NodeIndex(j);
@@ -347,6 +412,27 @@
 			newelements[i*this->elementswidth+j]=(int)sid+1;//C to Matlab indexing
 		}
-	}
- 
+		/*Verify the Jacobian determinant. If detJ<0, swap the 2 first postions:
+		  a -> b
+		  b -> a */
+		double detJ,xa,xb,xc,ya,yb,yc;
+		int a,b,c;
+
+		a=newelements[i*this->elementswidth+0]-1;
+		b=newelements[i*this->elementswidth+1]-1;
+		c=newelements[i*this->elementswidth+2]-1;
+
+		xa=newmeshX[a]; ya=newmeshY[a];
+		xb=newmeshX[b]; yb=newmeshY[b];
+		xc=newmeshX[c]; yc=newmeshY[c];
+
+		detJ=(xb-xa)*(yc-ya)-(xc-xa)*(yb-ya);
+	
+		/*verify and swap, if necessary*/
+		if(detJ<0) {
+			newelements[i*this->elementswidth+0]=b+1;//a->b
+			newelements[i*this->elementswidth+1]=a+1;//b->a
+		}
+	}
+
 	/*Setting outputs*/
 	nvertices	= nconformvertices;
@@ -358,5 +444,4 @@
 	/*Cleanup*/
 	xDelete<long>(vertex_index2sid);
-
 }
 /*}}}*/
