"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Kernel/GeoNormal.cpp" between
getdp-3.4.0-source.tgz and getdp-3.5.0-source.tgz

About: GetDP is a general finite element solver using mixed elements to discretize de Rham-type complexes in one, two and three dimensions.

GeoNormal.cpp  (getdp-3.4.0-source.tgz)	:	GeoNormal.cpp  (getdp-3.5.0-source.tgz)
// GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege		// GetDP - Copyright (C) 1997-2022 P. Dular and C. Geuzaine, University of Liege
//		//
// See the LICENSE.txt file for license information. Please report all		// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/getdp/getdp/issues.		// issues on https://gitlab.onelab.info/getdp/getdp/issues.
#include <string.h>		#include <string.h>
#include <math.h>		#include <math.h>
#include "ProData.h"		#include "ProData.h"
#include "GeoData.h"		#include "GeoData.h"
#include "MallocUtils.h"		#include "MallocUtils.h"
#include "Message.h"		#include "Message.h"
#define SQU(a) ((a)*(a))		#define SQU(a) ((a) * (a))
extern struct Problem Problem_S ;		extern struct Problem Problem_S;
extern struct GeoData * CurrentGeoData ;		extern struct GeoData *CurrentGeoData;
int fcmp_NXE(const void * a, const void * b)		int fcmp_NXE(const void *a, const void *b)
{		{
return		return ((struct Entity2XEntity1 *)a)->Num -
((struct Entity2XEntity1 *)a)->Num -		((struct Entity2XEntity1 *)b)->Num;
((struct Entity2XEntity1 *)b)->Num ;
}		}
int fcmp_EXVector(const void * a, const void * b)		int fcmp_EXVector(const void *a, const void *b)
{		{
return		return ((struct EntityXVector *)a)->Num - ((struct EntityXVector *)b)->Num;
((struct EntityXVector *)a)->Num -
((struct EntityXVector *)b)->Num ;
}		}
/*		/*
C'est une maniere un peu naive de creer cette BD. Mais elle a		C'est une maniere un peu naive de creer cette BD. Mais elle a
l'avantage de permettre une allocation simple (et minimum).		l'avantage de permettre une allocation simple (et minimum).
L'autre possibilite (boucler sur les elemnts) est plus rapide, mais		L'autre possibilite (boucler sur les elemnts) est plus rapide, mais
je ne vois pas bien comment obtenir un cout memoire minimum simplement,		je ne vois pas bien comment obtenir un cout memoire minimum simplement,
sans faire de nombreux realloc.		sans faire de nombreux realloc.
*/		*/
#define MAX_NBR_NXE_INCIDENCE 20		#define MAX_NBR_NXE_INCIDENCE 20
static int RegionIndexForNXE = -1;		static int RegionIndexForNXE = -1;
void Geo_CreateNodesXElements(int NumNode, int InIndex,		void Geo_CreateNodesXElements(int NumNode, int InIndex, int *NbrElements,
int *NbrElements, int **NumElements)		int **NumElements)
{		{
struct Entity2XEntity1 NXE, *NXE_P ;		struct Entity2XEntity1 NXE, *NXE_P;
struct Geo_Element *GeoElement ;		struct Geo_Element *GeoElement;
struct Group *Group_P ;		struct Group *Group_P;
int i, j, tmp[MAX_NBR_NXE_INCIDENCE] ;		int i, j, tmp[MAX_NBR_NXE_INCIDENCE];
Group_P = (struct Group*)List_Pointer(Problem_S.Group, InIndex);		Group_P = (struct Group *)List_Pointer(Problem_S.Group, InIndex);
if(InIndex != RegionIndexForNXE){		if(InIndex != RegionIndexForNXE) {
RegionIndexForNXE = InIndex ;		RegionIndexForNXE = InIndex;
Message::Info(" Generate NodesXElements information for Region '%s'", Group		Message::Info(" Generate NodesXElements information for Region '%s'",
_P->Name);		Group_P->Name);
if(CurrentGeoData->NodesXElements)		if(CurrentGeoData->NodesXElements)
Tree_Delete(CurrentGeoData->NodesXElements);		Tree_Delete(CurrentGeoData->NodesXElements);
CurrentGeoData->NodesXElements =		CurrentGeoData->NodesXElements =
Tree_Create(sizeof(struct Entity2XEntity1), fcmp_NXE) ;		Tree_Create(sizeof(struct Entity2XEntity1), fcmp_NXE);
}		}
NXE.Num = NumNode ;		NXE.Num = NumNode;
if((NXE_P = (struct Entity2XEntity1*)		if((NXE_P = (struct Entity2XEntity1 *)Tree_PQuery(
Tree_PQuery(CurrentGeoData->NodesXElements, &NXE))) {		CurrentGeoData->NodesXElements, &NXE))) {
*NbrElements = NXE_P->NbrEntities ;		*NbrElements = NXE_P->NbrEntities;
*NumElements = NXE_P->NumEntities ;		*NumElements = NXE_P->NumEntities;
}		}
else{		else {
NXE.NbrEntities = 0 ;		NXE.NbrEntities = 0;
for (i = 0 ; i < Geo_GetNbrGeoElements(); i++) {		for(i = 0; i < Geo_GetNbrGeoElements(); i++) {
GeoElement = Geo_GetGeoElement(i) ;		GeoElement = Geo_GetGeoElement(i);
if (List_Search(Group_P->InitialList, &GeoElement->Region, fcmp_int)){		if(List_Search(Group_P->InitialList, &GeoElement->Region, fcmp_int)) {
for(j=0 ; j<GeoElement->NbrNodes ; j++){		for(j = 0; j < GeoElement->NbrNodes; j++) {
if(GeoElement->NumNodes[j] == NumNode){		if(GeoElement->NumNodes[j] == NumNode) {
		/* printf("Adding elm %d to node %d\n", GeoElement->Num, NumNode);
/* printf("Adding elm %d to node %d\n", GeoElement->Num, NumNode); */		*/
/* this is to have orientation of elements adjacent to the node		/* this is to have orientation of elements adjacent to the node
Only valid for line elemnts !!!!!!! */		Only valid for line elemnts !!!!!!! */
tmp[NXE.NbrEntities] = ((!j)?-1:1) * GeoElement->Num ;		tmp[NXE.NbrEntities] = ((!j) ? -1 : 1) * GeoElement->Num;
NXE.NbrEntities++ ;		NXE.NbrEntities++;
}		}
}		}
}		}
}		}
NXE.NumEntities = (int*)Malloc(NXE.NbrEntities * sizeof(int)) ;		NXE.NumEntities = (int *)Malloc(NXE.NbrEntities * sizeof(int));
memcpy(NXE.NumEntities, tmp, NXE.NbrEntities * sizeof(int));		memcpy(NXE.NumEntities, tmp, NXE.NbrEntities * sizeof(int));
Tree_Add(CurrentGeoData->NodesXElements, &NXE);		Tree_Add(CurrentGeoData->NodesXElements, &NXE);
*NbrElements = NXE.NbrEntities ;		*NbrElements = NXE.NbrEntities;
*NumElements = NXE.NumEntities ;		*NumElements = NXE.NumEntities;
}		}
}		}
void Geo_CreateNormal(int Type, double *x, double *y, double *z, double *N)		void Geo_CreateNormal(int Type, double *x, double *y, double *z, double *N)
{		{
double x1x0, x2x0, y1y0, y2y0, z1z0, z2z0 ;		double x1x0, x2x0, y1y0, y2y0, z1z0, z2z0;
double nx, ny, nz, norm ;		double nx, ny, nz, norm;
switch (Type) {
case LINE :		switch(Type) {
case LINE_2 :		case LINE:
nx = y[1] - y[0] ;		case LINE_2:
ny = x[0] - x[1] ;		nx = y[1] - y[0];
norm = sqrt(SQU(nx)+SQU(ny)) ;		ny = x[0] - x[1];
N[0] = nx / norm ;		norm = sqrt(SQU(nx) + SQU(ny));
N[1] = ny / norm ;		N[0] = nx / norm;
N[2] = 0. ;		N[1] = ny / norm;
break ;		N[2] = 0.;
		break;
case TRIANGLE :
case TRIANGLE_2 :		case TRIANGLE:
case QUADRANGLE :		case TRIANGLE_2:
case QUADRANGLE_2 :		case QUADRANGLE:
x1x0 = x[1] - x[0] ;		case QUADRANGLE_2:
y1y0 = y[1] - y[0] ;		x1x0 = x[1] - x[0];
z1z0 = z[1] - z[0] ;		y1y0 = y[1] - y[0];
x2x0 = x[2] - x[0] ;		z1z0 = z[1] - z[0];
y2y0 = y[2] - y[0] ;		x2x0 = x[2] - x[0];
z2z0 = z[2] - z[0] ;		y2y0 = y[2] - y[0];
nx = y1y0 * z2z0 - z1z0 * y2y0 ;		z2z0 = z[2] - z[0];
ny = z1z0 * x2x0 - x1x0 * z2z0 ;		nx = y1y0 * z2z0 - z1z0 * y2y0;
nz = x1x0 * y2y0 - y1y0 * x2x0 ;		ny = z1z0 * x2x0 - x1x0 * z2z0;
norm = sqrt(SQU(nx)+SQU(ny)+SQU(nz)) ;		nz = x1x0 * y2y0 - y1y0 * x2x0;
N[0] = nx/norm ;		norm = sqrt(SQU(nx) + SQU(ny) + SQU(nz));
N[1] = ny/norm ;		N[0] = nx / norm;
N[2] = nz/norm ;		N[1] = ny / norm;
break ;		N[2] = nz / norm;
		break;
default :
Message::Error("Normal computation not done (yet) for Element Type %d", Type		default:
);		Message::Error("Normal computation not done (yet) for Element Type %d",
		Type);
}		}
}		}
void Geo_CreateNormalOfElement(struct Geo_Element *GeoElement, double *Normal)		void Geo_CreateNormalOfElement(struct Geo_Element *GeoElement, double *Normal)
{		{
struct EntityXVector EXV, *EXV_P ;		struct EntityXVector EXV, *EXV_P;
double x [NBR_MAX_NODES_IN_ELEMENT] ;		double x[NBR_MAX_NODES_IN_ELEMENT];
double y [NBR_MAX_NODES_IN_ELEMENT] ;		double y[NBR_MAX_NODES_IN_ELEMENT];
double z [NBR_MAX_NODES_IN_ELEMENT] ;		double z[NBR_MAX_NODES_IN_ELEMENT];
EXV.Num = GeoElement->Num ;		EXV.Num = GeoElement->Num;
if((EXV_P = (struct EntityXVector*)Tree_PQuery(CurrentGeoData->Normals, &EXV))		if((EXV_P =
) {		(struct EntityXVector *)Tree_PQuery(CurrentGeoData->Normals, &EXV))) {
memcpy(Normal, EXV_P->Vector, 3*sizeof(double));		memcpy(Normal, EXV_P->Vector, 3 * sizeof(double));
}		}
else{		else {
Geo_GetNodesCoordinates(GeoElement->NbrNodes, GeoElement->NumNodes, x, y, z)		Geo_GetNodesCoordinates(GeoElement->NbrNodes, GeoElement->NumNodes, x, y,
;		z);
Geo_CreateNormal(GeoElement->Type, x, y, z, Normal);		Geo_CreateNormal(GeoElement->Type, x, y, z, Normal);
memcpy(EXV.Vector, Normal, 3*sizeof(double));		memcpy(EXV.Vector, Normal, 3 * sizeof(double));
Tree_Add(CurrentGeoData->Normals, &EXV);		Tree_Add(CurrentGeoData->Normals, &EXV);
}		}
}		}
End of changes. 18 change blocks.
103 lines changed or deleted		99 lines changed or added

---