"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Functions/F_Gmsh.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.

[ To the main getdp source changes report ]

F_Gmsh.cpp  (getdp-3.4.0-source.tgz):F_Gmsh.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 "GetDPConfig.h" #include "GetDPConfig.h"
#include "ProData.h" #include "ProData.h"
#include "F.h" #include "F.h"
#include "Message.h" #include "Message.h"
extern struct CurrentData Current ; extern struct CurrentData Current;
#if defined(HAVE_GMSH) #if defined(HAVE_GMSH)
#include <gmsh/GmshGlobal.h> #include <gmsh/GmshGlobal.h>
#include <gmsh/PView.h> #include <gmsh/PView.h>
#include <gmsh/PViewData.h> #include <gmsh/PViewData.h>
#include <gmsh/Numeric.h> #include <gmsh/Numeric.h>
void F_Field(F_ARG) void F_Field(F_ARG)
{ {
if(A->Type != VECTOR){ if(A->Type != VECTOR) {
Message::Error("Field[] expects XYZ coordinates as argument"); Message::Error("Field[] expects XYZ coordinates as argument");
return; return;
} }
if(PView::list.empty()){ if(PView::list.empty()) {
Message::Error("No views available to interpolate from"); Message::Error("No views available to interpolate from");
return; return;
} }
double x = A->Val[0]; double x = A->Val[0];
double y = A->Val[1]; double y = A->Val[1];
double z = A->Val[2]; double z = A->Val[2];
if(Fct->NbrArguments > 1){ if(Fct->NbrArguments > 1) {
Message::Error("Time and additional arguments are not supported in Field: " Message::Error("Time and additional arguments are not supported in Field: "
"use {Scalar,Vector,Tensor}Field instead"); "use {Scalar,Vector,Tensor}Field instead");
return; return;
} }
for (int k = 0; k < Current.NbrHar; k++) for(int k = 0; k < Current.NbrHar; k++)
for (int j = 0; j < 9; j++) for(int j = 0; j < 9; j++) V->Val[MAX_DIM * k + j] = 0.;
V->Val[MAX_DIM * k + j] = 0. ;
V->Type = SCALAR; V->Type = SCALAR;
std::vector<int> iview; std::vector<int> iview;
if(!Fct->NbrParameters){ if(!Fct->NbrParameters) {
// use last view by default // use last view by default
iview.push_back(PView::list.back()->getTag()); iview.push_back(PView::list.back()->getTag());
} }
else{ else {
for(int i = 0; i < Fct->NbrParameters; i++) for(int i = 0; i < Fct->NbrParameters; i++) iview.push_back(Fct->Para[i]);
iview.push_back(Fct->Para[i]);
} }
double N = 0.; double N = 0.;
// add the values from all specified views // add the values from all specified views
for(unsigned int i = 0; i < iview.size(); i++){ for(unsigned int i = 0; i < iview.size(); i++) {
PView *v = PView::getViewByTag(iview[i]); PView *v = PView::getViewByTag(iview[i]);
if(!v){ if(!v) {
Message::Error("View with tag %d does not exist", iview[i]); Message::Error("View with tag %d does not exist", iview[i]);
return; return;
} }
PViewData *data = v->getData(); PViewData *data = v->getData();
std::vector<double> val(9 * data->getNumTimeSteps()); std::vector<double> val(9 * data->getNumTimeSteps());
if(data->searchScalar(x, y, z, &val[0])){ if(data->searchScalar(x, y, z, &val[0])) {
V->Val[0] += val[0]; V->Val[0] += val[0];
if(Current.NbrHar == 2 && data->getNumTimeSteps() > 1) if(Current.NbrHar == 2 && data->getNumTimeSteps() > 1)
V->Val[MAX_DIM] += val[1]; V->Val[MAX_DIM] += val[1];
V->Type = SCALAR; V->Type = SCALAR;
N += 1.; N += 1.;
} }
else if(data->searchVector(x, y, z, &val[0])){ else if(data->searchVector(x, y, z, &val[0])) {
for(int j = 0; j < 3; j++) for(int j = 0; j < 3; j++) V->Val[j] += val[j];
V->Val[j] += val[j]; if(Current.NbrHar == 2 && data->getNumTimeSteps() > 1) {
if(Current.NbrHar == 2 && data->getNumTimeSteps() > 1){ for(int j = 0; j < 3; j++) V->Val[MAX_DIM + j] += val[3 + j];
for(int j = 0; j < 3; j++)
V->Val[MAX_DIM + j] += val[3 + j];
} }
V->Type = VECTOR; V->Type = VECTOR;
N += 1.; N += 1.;
} }
else if(data->searchTensor(x, y, z, &val[0])){ else if(data->searchTensor(x, y, z, &val[0])) {
for(int j = 0; j < 9; j++) for(int j = 0; j < 9; j++) V->Val[j] += val[j];
V->Val[j] += val[j]; if(Current.NbrHar == 2 && data->getNumTimeSteps() > 1) {
if(Current.NbrHar == 2 && data->getNumTimeSteps() > 1){ for(int j = 0; j < 9; j++) V->Val[MAX_DIM + j] += val[9 + j];
for(int j = 0; j < 9; j++)
V->Val[MAX_DIM + j] += val[9 + j];
} }
V->Type = TENSOR; V->Type = TENSOR;
N += 1.; N += 1.;
} }
else{ else {
Message::Error("Did not find data at point (%g,%g,%g) in View with tag %d" Message::Error(
, "Did not find data at point (%g,%g,%g) in View with tag %d", x, y, z,
x, y, z, iview[i]); iview[i]);
} }
} }
if(N > 1.){ if(N > 1.) {
Message::Debug("Averaging data %g times on vertex (%g,%g,%g)", N, x, y, z); Message::Debug("Averaging data %g times on vertex (%g,%g,%g)", N, x, y, z);
for (int k = 0; k < Current.NbrHar; k++) for(int k = 0; k < Current.NbrHar; k++)
for (int j = 0; j < 9; j++) for(int j = 0; j < 9; j++) V->Val[MAX_DIM * k + j] = 0.;
V->Val[MAX_DIM * k + j] = 0. ;
} }
} }
static void F_X_Field(F_ARG, int type, bool complex, bool grad=false) static void F_X_Field(F_ARG, int type, bool complex, bool grad = false)
{ {
if(A->Type != VECTOR){ if(A->Type != VECTOR) {
Message::Error("Field[] expects XYZ coordinates as argument"); Message::Error("Field[] expects XYZ coordinates as argument");
return; return;
} }
if(PView::list.empty()){ if(PView::list.empty()) {
Message::Error("No views available to interpolate from"); Message::Error("No views available to interpolate from");
return; return;
} }
double x = A->Val[0]; double x = A->Val[0];
double y = A->Val[1]; double y = A->Val[1];
double z = A->Val[2]; double z = A->Val[2];
int numComp = (type == SCALAR) ? (grad ? 3 : 1) : int numComp = (type == SCALAR) ? (grad ? 3 : 1) :
(type == VECTOR) ? (grad ? 9 : 3) : 9; // TODO: grad of tensor (type == VECTOR) ? (grad ? 9 : 3) :
9; // TODO: grad of tensor
int NbrArg = Fct->NbrArguments ; int NbrArg = Fct->NbrArguments;
int TimeStep = 0, MatchElement = 0; int TimeStep = 0, MatchElement = 0;
if(NbrArg >= 2){ if(NbrArg >= 2) {
if((A+1)->Type != SCALAR){ if((A + 1)->Type != SCALAR) {
Message::Error("Expected scalar second argument (time step)"); Message::Error("Expected scalar second argument (time step)");
return; return;
} }
TimeStep = (int)(A+1)->Val[0]; TimeStep = (int)(A + 1)->Val[0];
} }
if(NbrArg >= 3){ if(NbrArg >= 3) {
if((A+2)->Type != SCALAR){ if((A + 2)->Type != SCALAR) {
Message::Error("Expected scalar second argument (element matching flag)"); Message::Error("Expected scalar second argument (element matching flag)");
return; return;
} }
MatchElement = (int)(A+2)->Val[0]; MatchElement = (int)(A + 2)->Val[0];
} }
// TODO: we could treat the third arguement as a tolerance (and call // TODO: we could treat the third arguement as a tolerance (and call
// searchScalarWithTol & friends) // searchScalarWithTol & friends)
// Complex{Scalar,Vector,Tensor}Field assume that the Gmsh view contains real // Complex{Scalar,Vector,Tensor}Field assume that the Gmsh view contains real
// and imaginary parts for each step // and imaginary parts for each step
if(complex) TimeStep *= 2; if(complex) TimeStep *= 2;
for (int k = 0; k < Current.NbrHar; k++) for(int k = 0; k < Current.NbrHar; k++)
for (int j = 0; j < numComp; j++) for(int j = 0; j < numComp; j++) V->Val[MAX_DIM * k + j] = 0.;
V->Val[MAX_DIM * k + j] = 0. ;
V->Type = (numComp == 1) ? SCALAR : (numComp == 3) ? VECTOR : TENSOR; V->Type = (numComp == 1) ? SCALAR : (numComp == 3) ? VECTOR : TENSOR;
std::vector<int> iview; std::vector<int> iview;
if(!Fct->NbrParameters){ if(!Fct->NbrParameters) {
// use last view by default // use last view by default
iview.push_back(PView::list.back()->getTag()); iview.push_back(PView::list.back()->getTag());
} }
else{ else {
for(int i = 0; i < Fct->NbrParameters; i++) for(int i = 0; i < Fct->NbrParameters; i++) iview.push_back(Fct->Para[i]);
iview.push_back(Fct->Para[i]);
} }
int qn = 0; int qn = 0;
double *qx = 0, *qy = 0, *qz = 0; double *qx = 0, *qy = 0, *qz = 0;
if(Current.Element){ if(Current.Element) {
qn = MatchElement ? Current.Element->GeoElement->NbrNodes : 0; qn = MatchElement ? Current.Element->GeoElement->NbrNodes : 0;
qx = Current.Element->x; qx = Current.Element->x;
qy = Current.Element->y; qy = Current.Element->y;
qz = Current.Element->z; qz = Current.Element->z;
} }
double N = 0.; double N = 0.;
// add the values from all specified views // add the values from all specified views
for(unsigned int i = 0; i < iview.size(); i++){ for(unsigned int i = 0; i < iview.size(); i++) {
PView *v = PView::getViewByTag(iview[i]); PView *v = PView::getViewByTag(iview[i]);
if(!v){ if(!v) {
Message::Error("View with tag %d does not exist", iview[i]); Message::Error("View with tag %d does not exist", iview[i]);
return; return;
} }
PViewData *data = v->getData(); PViewData *data = v->getData();
if(TimeStep < 0 || TimeStep >= data->getNumTimeSteps()){ if(TimeStep < 0 || TimeStep >= data->getNumTimeSteps()) {
Message::Error("Invalid step %d in View with tag %d", TimeStep, iview[i]); Message::Error("Invalid step %d in View with tag %d", TimeStep, iview[i]);
continue; continue;
} }
std::vector<double> val(numComp * data->getNumTimeSteps()); std::vector<double> val(numComp * data->getNumTimeSteps());
bool found = false; bool found = false;
switch(type){ switch(type) {
case SCALAR : case SCALAR:
if(data->searchScalar(x, y, z, &val[0], -1, 0, qn, qx, qy, qz, grad)) if(data->searchScalar(x, y, z, &val[0], -1, 0, qn, qx, qy, qz, grad))
found = true; found = true;
break; break;
case VECTOR : case VECTOR:
if(data->searchVector(x, y, z, &val[0], -1, 0, qn, qx, qy, qz, grad)) if(data->searchVector(x, y, z, &val[0], -1, 0, qn, qx, qy, qz, grad))
found = true; found = true;
break; break;
case TENSOR : case TENSOR:
// TODO: grad of tensor not allowed yet - not sure how we should return // TODO: grad of tensor not allowed yet - not sure how we should return
// the values; provide 3 routines that return 3 tensors, or add argumemt // the values; provide 3 routines that return 3 tensors, or add argumemt
// to select what to return? // to select what to return?
if(data->searchTensor(x, y, z, &val[0], -1, 0, qn, qx, qy, qz, false)) if(data->searchTensor(x, y, z, &val[0], -1, 0, qn, qx, qy, qz, false))
found = true; found = true;
break; break;
} }
if(found){ if(found) {
for(int j = 0; j < numComp; j++) for(int j = 0; j < numComp; j++) V->Val[j] += val[numComp * TimeStep + j];
V->Val[j] += val[numComp * TimeStep + j]; if(complex && Current.NbrHar == 2 &&
if(complex && Current.NbrHar == 2 && data->getNumTimeSteps() > TimeStep + data->getNumTimeSteps() > TimeStep + 1)
1)
for(int j = 0; j < numComp; j++) for(int j = 0; j < numComp; j++)
V->Val[MAX_DIM + j] += val[numComp * (TimeStep + 1) + j]; V->Val[MAX_DIM + j] += val[numComp * (TimeStep + 1) + j];
N += 1.; N += 1.;
} }
} }
if(N > 1.){ if(N > 1.) {
Message::Debug("Averaging data %g times on vertex (%g,%g,%g)", N, x, y, z); Message::Debug("Averaging data %g times on vertex (%g,%g,%g)", N, x, y, z);
for (int k = 0; k < Current.NbrHar; k++) for(int k = 0; k < Current.NbrHar; k++)
for (int j = 0; j < numComp; j++) for(int j = 0; j < numComp; j++) V->Val[MAX_DIM * k + j] /= N;
V->Val[MAX_DIM * k + j] /= N ;
} }
} }
void F_Distance(F_ARG) void F_Distance(F_ARG)
{ {
if(A->Type != VECTOR){ if(A->Type != VECTOR) {
Message::Error("Distance[] expects XYZ coordinates as argument"); Message::Error("Distance[] expects XYZ coordinates as argument");
return; return;
} }
SPoint3 p(A->Val[0], A->Val[1], A->Val[2]); SPoint3 p(A->Val[0], A->Val[1], A->Val[2]);
if(Fct->NbrParameters != 1){ if(Fct->NbrParameters != 1) {
Message::Error("Distance[] requires one view tag as parameter"); Message::Error("Distance[] requires one view tag as parameter");
return; return;
} }
PView *v = PView::getViewByTag(Fct->Para[0]); PView *v = PView::getViewByTag(Fct->Para[0]);
if(!v){ if(!v) {
Message::Error("View with tag %d does not exist", Fct->Para[0]); Message::Error("View with tag %d does not exist", Fct->Para[0]);
return; return;
} }
PViewData *data = v->getData(); PViewData *data = v->getData();
if(data->getNumTimeSteps() != 1 || !data->hasTimeStep(0)) { if(data->getNumTimeSteps() != 1 || !data->hasTimeStep(0)) {
Message::Error("View used in Distance[] should have a single step"); Message::Error("View used in Distance[] should have a single step");
return; return;
} }
skipping to change at line 265 skipping to change at line 256
if(data->skipEntity(0, ent)) continue; if(data->skipEntity(0, ent)) continue;
for(int ele = 0; ele < data->getNumElements(0, ent); ele++) { for(int ele = 0; ele < data->getNumElements(0, ent); ele++) {
if(data->skipElement(0, ent, ele)) continue; if(data->skipElement(0, ent, ele)) continue;
int numComp = data->getNumComponents(0, ent, ele); int numComp = data->getNumComponents(0, ent, ele);
if(numComp != 3) { if(numComp != 3) {
Message::Error("Distance[] expects a vector-valued view"); Message::Error("Distance[] expects a vector-valued view");
return; return;
} }
int numNodes = data->getNumNodes(0, ent, ele); int numNodes = data->getNumNodes(0, ent, ele);
if(numNodes != 2 && numNodes != 3) { if(numNodes != 2 && numNodes != 3) {
Message::Error("Distance[] can only currently handle lines and triangles Message::Error(
"); "Distance[] can only currently handle lines and triangles");
return; return;
} }
double nx[3], ny[3], nz[3], dx[3], dy[3], dz[3]; double nx[3], ny[3], nz[3], dx[3], dy[3], dz[3];
for(int nod = 0; nod < numNodes; nod++) { for(int nod = 0; nod < numNodes; nod++) {
data->getValue(0, ent, ele, nod, 0, dx[nod]); data->getValue(0, ent, ele, nod, 0, dx[nod]);
data->getValue(0, ent, ele, nod, 1, dy[nod]); data->getValue(0, ent, ele, nod, 1, dy[nod]);
data->getValue(0, ent, ele, nod, 2, dz[nod]); data->getValue(0, ent, ele, nod, 2, dz[nod]);
data->getNode(0, ent, ele, nod, nx[nod], ny[nod], nz[nod]); data->getNode(0, ent, ele, nod, nx[nod], ny[nod], nz[nod]);
} }
double d; double d;
skipping to change at line 291 skipping to change at line 283
signedDistancePointLine(p1, p2, p, d, closePt); signedDistancePointLine(p1, p2, p, d, closePt);
} }
else if(numNodes == 3) { else if(numNodes == 3) {
SPoint3 p3(nx[2] + dx[2], ny[2] + dy[2], nz[2] + dz[2]); SPoint3 p3(nx[2] + dx[2], ny[2] + dy[2], nz[2] + dz[2]);
signedDistancePointTriangle(p1, p2, p3, p, d, closePt); signedDistancePointTriangle(p1, p2, p3, p, d, closePt);
} }
dist = std::min(dist, d); dist = std::min(dist, d);
} }
} }
V->Type = SCALAR ; V->Type = SCALAR;
V->Val[0] = dist; V->Val[0] = dist;
V->Val[MAX_DIM] = 0.; V->Val[MAX_DIM] = 0.;
for (int k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) { for(int k = 2; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar); k += 2) {
V->Val[MAX_DIM* k] = V->Val[0] ; V->Val[MAX_DIM * k] = V->Val[0];
V->Val[MAX_DIM* (k+1)] = 0. ; V->Val[MAX_DIM * (k + 1)] = 0.;
} }
} }
#else #else
void F_Field(F_ARG) void F_Field(F_ARG)
{ {
Message::Error("You need to compile GetDP with Gmsh support to use Field[]"); Message::Error("You need to compile GetDP with Gmsh support to use Field[]");
V->Val[0] = 0. ; V->Val[0] = 0.;
V->Type = SCALAR ; V->Type = SCALAR;
} }
static void F_X_Field(F_ARG, int type, bool complex, bool grad=false) static void F_X_Field(F_ARG, int type, bool complex, bool grad = false)
{ {
Message::Error("You need to compile GetDP with Gmsh support to use Field[]"); Message::Error("You need to compile GetDP with Gmsh support to use Field[]");
V->Val[0] = 0. ; V->Val[0] = 0.;
V->Type = SCALAR ; V->Type = SCALAR;
} }
void F_Distance(F_ARG) void F_Distance(F_ARG)
{ {
Message::Error("You need to compile GetDP with Gmsh support to use Distance[]" Message::Error(
); "You need to compile GetDP with Gmsh support to use Distance[]");
V->Val[0] = 0. ; V->Val[0] = 0.;
V->Type = SCALAR ; V->Type = SCALAR;
} }
#endif #endif
void F_ScalarField(F_ARG){ F_X_Field(Fct, A, V, SCALAR, false); } void F_ScalarField(F_ARG) { F_X_Field(Fct, A, V, SCALAR, false); }
void F_VectorField(F_ARG){ F_X_Field(Fct, A, V, VECTOR, false); } void F_VectorField(F_ARG) { F_X_Field(Fct, A, V, VECTOR, false); }
void F_TensorField(F_ARG){ F_X_Field(Fct, A, V, TENSOR, false); } void F_TensorField(F_ARG) { F_X_Field(Fct, A, V, TENSOR, false); }
void F_ComplexScalarField(F_ARG){ F_X_Field(Fct, A, V, SCALAR, true); } void F_ComplexScalarField(F_ARG) { F_X_Field(Fct, A, V, SCALAR, true); }
void F_ComplexVectorField(F_ARG){ F_X_Field(Fct, A, V, VECTOR, true); } void F_ComplexVectorField(F_ARG) { F_X_Field(Fct, A, V, VECTOR, true); }
void F_ComplexTensorField(F_ARG){ F_X_Field(Fct, A, V, TENSOR, true); } void F_ComplexTensorField(F_ARG) { F_X_Field(Fct, A, V, TENSOR, true); }
void F_GradScalarField(F_ARG){ F_X_Field(Fct, A, V, SCALAR, false, true); } void F_GradScalarField(F_ARG) { F_X_Field(Fct, A, V, SCALAR, false, true); }
void F_GradVectorField(F_ARG){ F_X_Field(Fct, A, V, VECTOR, false, true); } void F_GradVectorField(F_ARG) { F_X_Field(Fct, A, V, VECTOR, false, true); }
void F_GradComplexScalarField(F_ARG){ F_X_Field(Fct, A, V, SCALAR, true, true); void F_GradComplexScalarField(F_ARG)
} {
void F_GradComplexVectorField(F_ARG){ F_X_Field(Fct, A, V, VECTOR, true, true); F_X_Field(Fct, A, V, SCALAR, true, true);
} }
void F_GradComplexVectorField(F_ARG)
{
F_X_Field(Fct, A, V, VECTOR, true, true);
}
 End of changes. 49 change blocks. 
89 lines changed or deleted 78 lines changed or added
To the C++ Programs section of the getdp source changes report or the top of this page
Mainly generated by pkgdiff using rfcdiff
Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)