From c74c22c96e87396f85d1473e4471633f280a6c6d Mon Sep 17 00:00:00 2001
From: Gauthier Becker <gauthierbecker@gmail.com>
Date: Thu, 29 Apr 2010 13:06:54 +0000
Subject: [PATCH] move contrib/cm3 to project/dgplate. Add FullDg formulation
Add Stress and InternalPoint (development)
---
CMakeLists.txt | 2 +-
Common/GmshMessage.cpp | 22 +-
Post/PViewDataGModelIO.cpp | 62 +-
Post/PViewDataListIO.cpp | 48 +-
Solver/dofManager.h | 5 +-
Solver/eigenSolver.cpp | 20 +-
Solver/eigenSolver.h | 4 +-
Solver/functionSpace.h | 3 +-
Solver/linearSystemCSR.cpp | 2 +-
Solver/linearSystemGMM.h | 9 +-
Solver/terms.h | 6 -
contrib/cm3/AssembleInterface.h | 54 --
contrib/cm3/C0DgPlateTerms.h | 650 -------------------
contrib/cm3/CMakeLists.txt | 8 -
contrib/cm3/DgC0PlateElementaryTerms.h | 387 ------------
contrib/cm3/DgC0PlateSolver.cpp | 386 ------------
contrib/cm3/DgC0PlateSolver.h | 96 ---
contrib/cm3/GModelWithInterface.cpp | 664 --------------------
contrib/cm3/GModelWithInterface.h | 59 --
contrib/cm3/LinearElasticShellHookeTensor.h | 75 ---
contrib/cm3/LocalBasis.cpp | 13 -
contrib/cm3/LocalBasis.h | 174 -----
contrib/cm3/MInterfaceElement.cpp | 13 -
contrib/cm3/MInterfaceElement.h | 167 -----
contrib/cm3/Tensor4dim2.h | 116 ----
contrib/cm3/mainDG.cpp | 39 --
26 files changed, 91 insertions(+), 2993 deletions(-)
delete mode 100644 contrib/cm3/AssembleInterface.h
delete mode 100644 contrib/cm3/C0DgPlateTerms.h
delete mode 100644 contrib/cm3/CMakeLists.txt
delete mode 100644 contrib/cm3/DgC0PlateElementaryTerms.h
delete mode 100644 contrib/cm3/DgC0PlateSolver.cpp
delete mode 100644 contrib/cm3/DgC0PlateSolver.h
delete mode 100644 contrib/cm3/GModelWithInterface.cpp
delete mode 100644 contrib/cm3/GModelWithInterface.h
delete mode 100644 contrib/cm3/LinearElasticShellHookeTensor.h
delete mode 100644 contrib/cm3/LocalBasis.cpp
delete mode 100644 contrib/cm3/LocalBasis.h
delete mode 100644 contrib/cm3/MInterfaceElement.cpp
delete mode 100644 contrib/cm3/MInterfaceElement.h
delete mode 100644 contrib/cm3/Tensor4dim2.h
delete mode 100644 contrib/cm3/mainDG.cpp
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 43c89bee2d..ca602e655a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -866,7 +866,7 @@ if(ENABLE_ARC)
include(contrib/arc/CMakeLists.txt)
endif(ENABLE_ARC)
if(ENABLE_CM3)
- include(contrib/cm3/CMakeLists.txt)
+ include(projects/dgplate/CMakeLists.txt)
endif(ENABLE_CM3)
find_program(BISON bison)
diff --git a/Common/GmshMessage.cpp b/Common/GmshMessage.cpp
index 0ba814037c..26a052d022 100644
--- a/Common/GmshMessage.cpp
+++ b/Common/GmshMessage.cpp
@@ -115,10 +115,10 @@ void Msg::Exit(int level)
// the persistent info to disk
if(FlGui::available() && !_commRank) {
if(CTX::instance()->sessionSave)
- PrintOptions(0, GMSH_SESSIONRC, 0, 0,
+ PrintOptions(0, GMSH_SESSIONRC, 0, 0,
(CTX::instance()->homeDir + CTX::instance()->sessionFileName).c_str());
if(CTX::instance()->optionsSave)
- PrintOptions(0, GMSH_OPTIONSRC, 1, 0,
+ PrintOptions(0, GMSH_OPTIONSRC, 1, 0,
(CTX::instance()->homeDir + CTX::instance()->optionsFileName).c_str());
}
#endif
@@ -199,7 +199,7 @@ void Msg::Error(const char *fmt, ...)
#endif
if(CTX::instance()->terminal){
- if(_commSize > 1)
+ if(_commSize > 1)
fprintf(stderr, "Error : [On processor %d] %s\n", _commRank, str);
else
fprintf(stderr, "Error : %s\n", str);
@@ -302,7 +302,7 @@ void Msg::Direct(int level, const char *fmt, ...)
FlGui::instance()->messages->show();
}
#endif
-
+
if(CTX::instance()->terminal){
fprintf(stdout, "%s\n", str);
fflush(stdout);
@@ -326,7 +326,7 @@ void Msg::StatusBar(int num, bool log, const char *fmt, ...)
#if defined(HAVE_FLTK)
if(FlGui::available()){
if(log) FlGui::instance()->check();
- if(!log || num != 2 || _verbosity > 3)
+ if(!log || num != 2 || _verbosity > 3)
FlGui::instance()->setStatus(str, num - 1);
if(log){
std::string tmp = std::string("Info : ") + str;
@@ -362,7 +362,7 @@ void Msg::Debug(const char *fmt, ...)
#endif
if(CTX::instance()->terminal){
- if(_commSize > 1)
+ if(_commSize > 1)
fprintf(stdout, "Debug : [On processor %d] %s\n", _commRank, str);
else
fprintf(stdout, "Debug : %s\n", str);
@@ -427,7 +427,7 @@ void Msg::PrintTimers()
{
// do a single stdio call!
std::string str;
- for(std::map<std::string, double>::iterator it = _timers.begin();
+ for(std::map<std::string, double>::iterator it = _timers.begin();
it != _timers.end(); it++){
if(it != _timers.begin()) str += ", ";
char tmp[256];
@@ -437,7 +437,7 @@ void Msg::PrintTimers()
if(!str.size()) return;
if(CTX::instance()->terminal){
- if(_commSize > 1)
+ if(_commSize > 1)
fprintf(stdout, "Timers : [On processor %d] %s\n", _commRank, str.c_str());
else
fprintf(stdout, "Timers : %s\n", str.c_str());
@@ -474,9 +474,9 @@ void Msg::PrintErrorCounter(const char *title)
#endif
if(CTX::instance()->terminal){
- fprintf(stderr, "%s\n%s\n%s\n%s\n%s\n%s\n", (prefix + line).c_str(),
+ fprintf(stderr, "%s\n%s\n%s\n%s\n%s\n%s\n", (prefix + line).c_str(),
(prefix + title).c_str(), (prefix + warn).c_str(),
- (prefix + err).c_str(), (prefix + help).c_str(),
+ (prefix + err).c_str(), (prefix + help).c_str(),
(prefix + line).c_str());
fflush(stderr);
}
@@ -509,7 +509,7 @@ double Msg::GetValue(const char *text, double defaultval)
return atof(str);
}
-int Msg::GetAnswer(const char *question, int defaultval, const char *zero,
+int Msg::GetAnswer(const char *question, int defaultval, const char *zero,
const char *one, const char *two)
{
// if a callback is given let's assume we don't want to be bothered
diff --git a/Post/PViewDataGModelIO.cpp b/Post/PViewDataGModelIO.cpp
index 97a266e439..3eb0a19aa8 100644
--- a/Post/PViewDataGModelIO.cpp
+++ b/Post/PViewDataGModelIO.cpp
@@ -18,7 +18,7 @@ bool PViewDataGModel::addData(GModel *model, std::map<int, std::vector<double> >
int numComp = 9;
if (numC < 0){
- for(std::map<int, std::vector<double> >::iterator it = data.begin();
+ for(std::map<int, std::vector<double> >::iterator it = data.begin();
it != data.end(); it++)
numComp = std::min(numComp, (int)it->second.size());
}
@@ -26,14 +26,14 @@ bool PViewDataGModel::addData(GModel *model, std::map<int, std::vector<double> >
while(step >= (int)_steps.size())
_steps.push_back(new stepData<double>(model, numComp));
-
+
_steps[step]->setTime(time);
-
- int numEnt = (_type == NodeData) ? model->getNumMeshVertices() :
+
+ int numEnt = (_type == NodeData) ? model->getNumMeshVertices() :
model->getNumMeshElements();
_steps[step]->resizeData(numEnt);
-
- for(std::map<int, std::vector<double> >::iterator it = data.begin();
+
+ for(std::map<int, std::vector<double> >::iterator it = data.begin();
it != data.end(); it++){
double *d = _steps[step]->getData(it->first, true);
for(int j = 0; j < numComp; j++)
@@ -45,15 +45,15 @@ bool PViewDataGModel::addData(GModel *model, std::map<int, std::vector<double> >
}
bool PViewDataGModel::readMSH(std::string fileName, int fileIndex, FILE *fp,
- bool binary, bool swap, int step, double time,
+ bool binary, bool swap, int step, double time,
int partition, int numComp, int numEnt)
{
- Msg::Info("Reading step %d (time %g) partition %d: %d records",
+ Msg::Info("Reading step %d (time %g) partition %d: %d records",
step, time, partition, numEnt);
while(step >= (int)_steps.size())
_steps.push_back(new stepData<double>(GModel::current(), numComp));
-
+
_steps[step]->setFileName(fileName);
_steps[step]->setFileIndex(fileIndex);
_steps[step]->setTime(time);
@@ -88,7 +88,7 @@ bool PViewDataGModel::readMSH(std::string fileName, int fileIndex, FILE *fp,
}
double *d = _steps[step]->getData(num, true, mult);
if(binary){
- if((int)fread(d, sizeof(double), numComp * mult, fp) != numComp * mult)
+ if((int)fread(d, sizeof(double), numComp * mult, fp) != numComp * mult)
return false;
if(swap) SwapBytes((char*)d, sizeof(double), numComp * mult);
}
@@ -96,7 +96,7 @@ bool PViewDataGModel::readMSH(std::string fileName, int fileIndex, FILE *fp,
for(int j = 0; j < numComp * mult; j++)
if(fscanf(fp, "%lf", &d[j]) != 1) return false;
}
- if(numEnt > 100000)
+ if(numEnt > 100000)
Msg::ProgressMeter(i + 1, numEnt, "Reading data");
}
@@ -209,7 +209,7 @@ bool PViewDataGModel::writeMSH(std::string fileName, bool binary)
}
}
}
-
+
fclose(fp);
return true;
}
@@ -230,7 +230,7 @@ bool PViewDataGModel::readMED(std::string fileName, int fileIndex)
Msg::Error("Unable to open file '%s'", fileName.c_str());
return false;
}
-
+
med_int numComp = MEDnChamp(fid, fileIndex + 1);
if(numComp <= 0){
Msg::Error("Could not get number of components for MED field");
@@ -241,7 +241,7 @@ bool PViewDataGModel::readMED(std::string fileName, int fileIndex)
std::vector<char> compName(numComp * MED_TAILLE_PNOM + 1);
std::vector<char> compUnit(numComp * MED_TAILLE_PNOM + 1);
med_type_champ type;
- if(MEDchampInfo(fid, fileIndex + 1, name, &type, &compName[0], &compUnit[0],
+ if(MEDchampInfo(fid, fileIndex + 1, name, &type, &compName[0], &compUnit[0],
numComp) < 0){
Msg::Error("Could not get MED field info");
return false;
@@ -250,7 +250,7 @@ bool PViewDataGModel::readMED(std::string fileName, int fileIndex)
Msg::Info("Reading %d-component field <<%s>>", numComp, name);
setName(name);
- int numCompMsh =
+ int numCompMsh =
(numComp <= 1) ? 1 : (numComp <= 3) ? 3 : (numComp <= 9) ? 9 : numComp;
if(numCompMsh > 9) Msg::Warning("More than 9 components in field");
@@ -259,13 +259,13 @@ bool PViewDataGModel::readMED(std::string fileName, int fileIndex)
// important: it should match the ordering of the MSH element types
// (when elements are saved without tags, this governs the order
// with which we implicitly index them in GModel::readMED)
- const med_entite_maillage entType[] =
+ const med_entite_maillage entType[] =
{MED_NOEUD, MED_MAILLE, MED_NOEUD_MAILLE};
- const med_geometrie_element eleType[] =
- {MED_NONE, MED_SEG2, MED_TRIA3, MED_QUAD4, MED_TETRA4, MED_HEXA8,
- MED_PENTA6, MED_PYRA5, MED_SEG3, MED_TRIA6, MED_TETRA10,
+ const med_geometrie_element eleType[] =
+ {MED_NONE, MED_SEG2, MED_TRIA3, MED_QUAD4, MED_TETRA4, MED_HEXA8,
+ MED_PENTA6, MED_PYRA5, MED_SEG3, MED_TRIA6, MED_TETRA10,
MED_POINT1, MED_QUAD8, MED_HEXA20, MED_PENTA15, MED_PYRA13};
- const int nodesPerEle[] =
+ const int nodesPerEle[] =
{0, 2, 3, 4, 4, 8, 6, 5, 3, 6, 10, 1, 8, 20, 15, 13};
med_int numSteps = 0;
@@ -279,7 +279,7 @@ bool PViewDataGModel::readMED(std::string fileName, int fileIndex)
}
if(!i && !j) break; // MED_NOEUD does not care about eleType
}
- }
+ }
if(numSteps < 1 || pairs.empty()){
Msg::Error("Nothing to import from MED file");
@@ -287,7 +287,7 @@ bool PViewDataGModel::readMED(std::string fileName, int fileIndex)
}
else{
med_entite_maillage ent = entType[pairs[0].first];
- _type = (ent == MED_NOEUD) ? NodeData : (ent == MED_MAILLE) ? ElementData :
+ _type = (ent == MED_NOEUD) ? NodeData : (ent == MED_MAILLE) ? ElementData :
ElementNodeData;
}
@@ -340,7 +340,7 @@ bool PViewDataGModel::readMED(std::string fileName, int fileIndex)
std::vector<double> val(numVal * numComp);
char locname[MED_TAILLE_NOM + 1], profileName[MED_TAILLE_NOM + 1];
if(MEDchampLire(fid, meshName, name, (unsigned char*)&val[0], MED_FULL_INTERLACE,
- MED_ALL, locname, profileName, MED_COMPACT, ent, ele,
+ MED_ALL, locname, profileName, MED_COMPACT, ent, ele,
numdt, numo) < 0){
Msg::Error("Could not read field values");
return false;
@@ -367,7 +367,7 @@ bool PViewDataGModel::readMED(std::string fileName, int fileIndex)
// internal reference element
for(int i = 0; i < (int)gscoo.size(); i++){
p.push_back(gscoo[i]);
- if(i % dim == dim - 1) for(int j = 0; j < 3 - dim; j++) p.push_back(0.);
+ if(i % dim == dim - 1) for(int j = 0; j < 3 - dim; j++) p.push_back(0.);
}
}
}
@@ -392,9 +392,9 @@ bool PViewDataGModel::readMED(std::string fileName, int fileIndex)
// get size of full array and tags (if any) of entities
bool nodal = (ent == MED_NOEUD);
- med_int numEnt = MEDnEntMaa(fid, meshName, nodal ? MED_COOR : MED_CONN,
- nodal ? MED_NOEUD : MED_MAILLE,
- nodal ? MED_NONE : ele,
+ med_int numEnt = MEDnEntMaa(fid, meshName, nodal ? MED_COOR : MED_CONN,
+ nodal ? MED_NOEUD : MED_MAILLE,
+ nodal ? MED_NONE : ele,
nodal ? (med_connectivite)0 : MED_NOD);
std::vector<med_int> tags(numEnt);
if(MEDnumLire(fid, meshName, &tags[0], numEnt, nodal ? MED_NOEUD : MED_MAILLE,
@@ -506,7 +506,7 @@ bool PViewDataGModel::writeMED(std::string fileName)
return false;
}
- med_int numNodes = MEDnEntMaa(fid, meshName, MED_COOR, MED_NOEUD,
+ med_int numNodes = MEDnEntMaa(fid, meshName, MED_COOR, MED_NOEUD,
MED_NONE, (med_connectivite)0);
if(numNodes <= 0){
Msg::Error("Could not get valid number of nodes in mesh");
@@ -525,7 +525,7 @@ bool PViewDataGModel::writeMED(std::string fileName)
for(unsigned int i = 0; i < profile.size(); i++)
for(int k = 0; k < numComp; k++)
val[i * numComp + k] = _steps[step]->getData(indices[i])[k];
- if(MEDchampEcr(fid, meshName, fieldName, (unsigned char*)&val[0],
+ if(MEDchampEcr(fid, meshName, fieldName, (unsigned char*)&val[0],
MED_FULL_INTERLACE, numNodes, (char*)MED_NOGAUSS, MED_ALL,
profileName, MED_COMPACT, MED_NOEUD, MED_NONE, (med_int)step,
(char*)"unknown", time, MED_NONOR) < 0) {
@@ -533,7 +533,7 @@ bool PViewDataGModel::writeMED(std::string fileName)
return false;
}
}
-
+
if(MEDfermer(fid) < 0){
Msg::Error("Unable to close file '%s'", (char*)fileName.c_str());
return false;
@@ -545,7 +545,7 @@ bool PViewDataGModel::writeMED(std::string fileName)
bool PViewDataGModel::readMED(std::string fileName, int fileIndex)
{
- Msg::Error("Gmsh must be compiled with MED support to read '%s'",
+ Msg::Error("Gmsh must be compiled with MED support to read '%s'",
fileName.c_str());
return false;
}
diff --git a/Post/PViewDataListIO.cpp b/Post/PViewDataListIO.cpp
index d32ac0ac77..9eba8db064 100644
--- a/Post/PViewDataListIO.cpp
+++ b/Post/PViewDataListIO.cpp
@@ -13,7 +13,7 @@
#include "Context.h"
#include "adaptiveData.h"
-static void dVecRead(std::vector<double> &v, int n, FILE *fp,
+static void dVecRead(std::vector<double> &v, int n, FILE *fp,
bool binary, int swap)
{
if(!n) return;
@@ -33,7 +33,7 @@ static void dVecRead(std::vector<double> &v, int n, FILE *fp,
}
}
-static void cVecRead(std::vector<char> &v, int n, FILE *fp,
+static void cVecRead(std::vector<char> &v, int n, FILE *fp,
bool binary, int swap, bool oldStyle)
{
if(!n) return;
@@ -115,7 +115,7 @@ bool PViewDataList::readPOS(FILE *fp, double version, bool binary)
else if(version == 1.1) {
Msg::Debug("Detected post-processing view format 1.1");
if(!fscanf(fp, "%s %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",
- name, &NbTimeStep, &NbSP, &NbVP, &NbTP, &NbSL, &NbVL, &NbTL,
+ name, &NbTimeStep, &NbSP, &NbVP, &NbTP, &NbSL, &NbVL, &NbTL,
&NbST, &NbVT, &NbTT, &NbSS, &NbVS, &NbTS, &NbT2, &t2l, &NbT3,
&t3l)){
Msg::Error("Read error");
@@ -127,7 +127,7 @@ bool PViewDataList::readPOS(FILE *fp, double version, bool binary)
if(!fscanf(fp, "%s %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d "
"%d %d %d %d %d %d %d %d %d %d %d %d %d\n",
name, &NbTimeStep, &NbSP, &NbVP, &NbTP, &NbSL, &NbVL, &NbTL,
- &NbST, &NbVT, &NbTT, &NbSQ, &NbVQ, &NbTQ, &NbSS, &NbVS, &NbTS,
+ &NbST, &NbVT, &NbTT, &NbSQ, &NbVQ, &NbTQ, &NbSS, &NbVS, &NbTS,
&NbSH, &NbVH, &NbTH, &NbSI, &NbVI, &NbTI, &NbSY, &NbVY, &NbTY,
&NbT2, &t2l, &NbT3, &t3l)){
Msg::Error("Read error");
@@ -140,10 +140,10 @@ bool PViewDataList::readPOS(FILE *fp, double version, bool binary)
"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d "
"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",
name, &NbTimeStep, &NbSP, &NbVP, &NbTP, &NbSL, &NbVL, &NbTL,
- &NbST, &NbVT, &NbTT, &NbSQ, &NbVQ, &NbTQ, &NbSS, &NbVS, &NbTS,
+ &NbST, &NbVT, &NbTT, &NbSQ, &NbVQ, &NbTQ, &NbSS, &NbVS, &NbTS,
&NbSH, &NbVH, &NbTH, &NbSI, &NbVI, &NbTI, &NbSY, &NbVY, &NbTY,
- &NbSL2, &NbVL2, &NbTL2, &NbST2, &NbVT2, &NbTT2, &NbSQ2, &NbVQ2,
- &NbTQ2, &NbSS2, &NbVS2, &NbTS2, &NbSH2, &NbVH2, &NbTH2, &NbSI2,
+ &NbSL2, &NbVL2, &NbTL2, &NbST2, &NbVT2, &NbTT2, &NbSQ2, &NbVQ2,
+ &NbTQ2, &NbSS2, &NbVS2, &NbTS2, &NbSH2, &NbVH2, &NbTH2, &NbSI2,
&NbVI2, &NbTI2, &NbSY2, &NbVY2, &NbTY2, &NbT2, &t2l, &NbT3, &t3l)){
Msg::Error("Read error");
return false;
@@ -157,7 +157,7 @@ bool PViewDataList::readPOS(FILE *fp, double version, bool binary)
for(int i = 0; i < (int)strlen(name); i++)
if(name[i] == '^')
name[i] = ' ';
-
+
int swap = 0;
if(binary) {
int testone;
@@ -170,7 +170,7 @@ bool PViewDataList::readPOS(FILE *fp, double version, bool binary)
swap = 1;
}
}
-
+
dVecRead(Time, NbTimeStep, fp, binary, swap);
dVecRead(SP, NbSP * (NbTimeStep * 1 + 3), fp, binary, swap);
dVecRead(VP, NbVP * (NbTimeStep * 3 + 3), fp, binary, swap);
@@ -240,12 +240,12 @@ bool PViewDataList::readPOS(FILE *fp, double version, bool binary)
if(NbSY2){ NbSY = NbSY2; setOrder2(TYPE_PYR); }
if(NbVY2){ NbVY = NbVY2; setOrder2(TYPE_PYR); }
if(NbTY2){ NbTY = NbTY2; setOrder2(TYPE_PYR); }
-
+
dVecRead(T2D, NbT2 * 4, fp, binary, swap);
cVecRead(T2C, t2l, fp, binary, swap, (version <= 1.2));
dVecRead(T3D, NbT3 * 5, fp, binary, swap);
cVecRead(T3C, t3l, fp, binary, swap, (version <= 1.2));
-
+
Msg::Debug("Read View '%s' (%d TimeSteps): "
"SP(%d/%d) VP(%d/%d) TP(%d/%d) "
"SL(%d/%d) VL(%d/%d) TL(%d/%d) "
@@ -254,8 +254,8 @@ bool PViewDataList::readPOS(FILE *fp, double version, bool binary)
"SS(%d/%d) VS(%d/%d) TS(%d/%d) "
"SH(%d/%d) VH(%d/%d) TH(%d/%d) "
"SI(%d/%d) VI(%d/%d) TI(%d/%d) "
- "SY(%d/%d) VY(%d/%d) TY(%d/%d) "
- "T2(%d/%d/%d) T3(%d/%d/%d) ",
+ "SY(%d/%d) VY(%d/%d) TY(%d/%d) "
+ "T2(%d/%d/%d) T3(%d/%d/%d) ",
name, NbTimeStep,
NbSP, SP.size(), NbVP, VP.size(), NbTP, TP.size(),
NbSL, SL.size(), NbVL, VL.size(), NbTL, TL.size(),
@@ -265,9 +265,9 @@ bool PViewDataList::readPOS(FILE *fp, double version, bool binary)
NbSH, SH.size(), NbVH, VH.size(), NbTH, TH.size(),
NbSI, SI.size(), NbVI, VI.size(), NbTI, TI.size(),
NbSY, SY.size(), NbVY, VY.size(), NbTY, TY.size(),
- NbT2, T2D.size(), T2C.size(),
+ NbT2, T2D.size(), T2C.size(),
NbT3, T3D.size(), T3C.size());
-
+
setName(name);
finalize();
return true;
@@ -309,7 +309,7 @@ static void writeElementPOS(FILE *fp, const char *str, int nbnod, int nb,
}
}
-static void writeTextPOS(FILE *fp, int nbc, int nb, std::vector<double> &TD,
+static void writeTextPOS(FILE *fp, int nbc, int nb, std::vector<double> &TD,
std::vector<char> &TC)
{
if(!nb || (nbc != 4 && nbc != 5)) return;
@@ -383,7 +383,7 @@ bool PViewDataList::writePOS(std::string fileName, bool binary, bool parsed, boo
}
}
dVecWrite(Time, fp, binary);
- dVecWrite(SP, fp, binary); dVecWrite(VP, fp, binary);
+ dVecWrite(SP, fp, binary); dVecWrite(VP, fp, binary);
dVecWrite(TP, fp, binary); dVecWrite(SL, fp, binary);
dVecWrite(VL, fp, binary); dVecWrite(TL, fp, binary);
dVecWrite(ST, fp, binary); dVecWrite(VT, fp, binary);
@@ -446,7 +446,7 @@ class pVertexLessThan{
}
};
-static void getNodeMSH(int nbelm, std::vector<double> &list,
+static void getNodeMSH(int nbelm, std::vector<double> &list,
int nbnod, int nbcomp, int nbstep,
std::set<pVertex, pVertexLessThan> *nodes, int *numelm)
{
@@ -462,8 +462,8 @@ static void getNodeMSH(int nbelm, std::vector<double> &list,
std::set<pVertex, pVertexLessThan>::iterator it = nodes->find(n);
if(it == nodes->end()){
n.Num = nodes->size() + 1;
- for(int ts = 0; ts < nbstep; ts++)
- for(int k = 0; k < nbcomp; k++)
+ for(int ts = 0; ts < nbstep; ts++)
+ for(int k = 0; k < nbcomp; k++)
n.Val.push_back(v[nbcomp * nbnod * ts + nbcomp * j + k]);
nodes->insert(n);
}
@@ -472,7 +472,7 @@ static void getNodeMSH(int nbelm, std::vector<double> &list,
}
}
-static void writeElementMSH(FILE *fp, int num, int nbnod, pVertex nod[8],
+static void writeElementMSH(FILE *fp, int num, int nbnod, pVertex nod[8],
int nbcomp, double *vals, int dim)
{
switch(dim){
@@ -486,7 +486,7 @@ static void writeElementMSH(FILE *fp, int num, int nbnod, pVertex nod[8],
if(nbnod == 3)
fprintf(fp, "%d 2 0 %d %d %d\n", num, nod[0].Num, nod[1].Num, nod[2].Num);
else
- fprintf(fp, "%d 3 0 %d %d %d %d\n", num, nod[0].Num, nod[1].Num,
+ fprintf(fp, "%d 3 0 %d %d %d %d\n", num, nod[0].Num, nod[1].Num,
nod[2].Num, nod[3].Num);
break;
case 3:
@@ -501,7 +501,7 @@ static void writeElementMSH(FILE *fp, int num, int nbnod, pVertex nod[8],
fprintf(fp, "%d 6 0 %d %d %d %d %d %d\n", num, nod[0].Num, nod[1].Num,
nod[2].Num, nod[3].Num, nod[4].Num, nod[5].Num);
else
- fprintf(fp, "%d 5 0 %d %d %d %d %d %d %d %d\n", num, nod[0].Num,
+ fprintf(fp, "%d 5 0 %d %d %d %d %d %d %d %d\n", num, nod[0].Num,
nod[1].Num, nod[2].Num, nod[3].Num, nod[4].Num, nod[5].Num,
nod[6].Num, nod[7].Num);
break;
@@ -509,7 +509,7 @@ static void writeElementMSH(FILE *fp, int num, int nbnod, pVertex nod[8],
}
static void writeElementsMSH(FILE *fp, int nbelm, std::vector<double> &list,
- int nbnod, int nbcomp, int dim,
+ int nbnod, int nbcomp, int dim,
std::set<pVertex, pVertexLessThan> *nodes,
int *numelm)
{
diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index dfeb3f4bf6..08b56f44af 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -17,7 +17,7 @@
#include <iostream>
class Dof{
- private:
+ protected:
// v(x) = \sum_f \sum_i v_{fi} s^f_i(x)
long int _entity; // "i": node, edge, group, etc.
int _type; // "f": basis function type index, etc.
@@ -34,6 +34,9 @@ class Dof{
i1 = t % 10000;
i2 = t / 10000;
}
+ //inline static int createEntityForFullDg(const int ele_num, const int ver_num){return ele_num*100000+ver_num;} // Problem if the number of Vertex exceed 1e5 OK for a moment
+ //inline static void getTwoIntFromEntityFullDg(int t, int &i1, int &i2){i1 = t % 100000 ; i2 = t / 100000;}
+
bool operator < (const Dof &other) const
{
if(_entity < other._entity) return true;
diff --git a/Solver/eigenSolver.cpp b/Solver/eigenSolver.cpp
index 4e58cbad6c..dac0349778 100644
--- a/Solver/eigenSolver.cpp
+++ b/Solver/eigenSolver.cpp
@@ -10,7 +10,7 @@
#include <slepceps.h>
-eigenSolver::eigenSolver(dofManager<double> *manager, std::string A,
+eigenSolver::eigenSolver(dofManager<double> *manager, std::string A,
std::string B, bool hermitian) : _A(0), _B(0),_hermitian(hermitian)
{
if(A.size()){
@@ -29,7 +29,7 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
if(!_A) return false;
Mat A = _A->getMatrix();
Mat B = _B ? _B->getMatrix() : PETSC_NULL;
-
+
_try(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
_try(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
PetscInt N, M;
@@ -40,7 +40,7 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
_try(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
_try(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
_try(MatGetSize(B, &N2, &M2));
- }
+ }
// generalized eigenvalue problem A x - \lambda B x = 0
EPS eps;
@@ -56,7 +56,7 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
_try(EPSSetTolerances(eps,1.e-6,20));//1.e-7 50
//_try(EPSSetType(eps, EPSKRYLOVSCHUR)); //default
_try(EPSSetType(eps,EPSARNOLDI));
- //_try(EPSSetType(eps,EPSARPACK));
+ //_try(EPSSetType(eps,EPSARPACK));
//_try(EPSSetType(eps,EPSPOWER));
// override these options at runtime, petsc-style
@@ -89,7 +89,7 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
Msg::Info("SLEPc solving...");
double t1 = Cpu();
_try(EPSSolve(eps));
-
+
// check convergence
int its;
_try(EPSGetIterationNumber(eps, &its));
@@ -105,7 +105,7 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
Msg::Error("SLEPc generic breakdown in method");
else if(reason == EPS_DIVERGED_NONSYMMETRIC)
Msg::Error("The operator is nonsymmetric");
-
+
// get number of converged approximate eigenpairs
PetscInt nconv;
_try(EPSGetConverged(eps, &nconv));
@@ -113,7 +113,7 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
// ignore additional eigenvalues if we get more than what we asked
if(nconv > nev) nconv = nev;
-
+
if (nconv > 0) {
Vec xr, xi;
_try(MatGetVecs(A, PETSC_NULL, &xr));
@@ -132,9 +132,9 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
PetscReal re = kr;
PetscReal im = ki;
#endif
- Msg::Info("EIG %03d %s%.16e %s%.16e %3.6e",
+ Msg::Info("EIG %03d %s%.16e %s%.16e %3.6e",
i, (re < 0) ? "" : " ", re, (im < 0) ? "" : " ", im, error);
-
+
// store eigenvalues and eigenvectors
_eigenValues.push_back(std::complex<double>(re, im));
PetscScalar *tmpr, *tmpi;
@@ -156,7 +156,7 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
_try(EPSDestroy(eps));
_try(SlepcFinalize());
-
+
if(reason == EPS_CONVERGED_TOL){
Msg::Info("SLEPc done");
return true;
diff --git a/Solver/eigenSolver.h b/Solver/eigenSolver.h
index 677e35a582..1224f39ced 100644
--- a/Solver/eigenSolver.h
+++ b/Solver/eigenSolver.h
@@ -25,7 +25,7 @@ class eigenSolver{
std::vector<std::vector<std::complex<double> > > _eigenVectors;
void _try(int ierr) const { CHKERRABORT(PETSC_COMM_WORLD, ierr); }
public:
- eigenSolver(dofManager<double> *manager, std::string A,
+ eigenSolver(dofManager<double> *manager, std::string A,
std::string B="", bool hermitian=false);
bool solve(int numEigenValues=0, std::string which="");
int getNumEigenValues(){ return _eigenValues.size(); }
@@ -39,7 +39,7 @@ class eigenSolver{
private:
std::vector<std::complex<double> > _dummy;
public:
- eigenSolver(dofManager<double> *manager, std::string A,
+ eigenSolver(dofManager<double> *manager, std::string A,
std::string B="", bool hermitian=false){}
bool solve(int numEigenValues=0, std::string which="")
{
diff --git a/Solver/functionSpace.h b/Solver/functionSpace.h
index 13b7d83f0a..53ea4afd43 100644
--- a/Solver/functionSpace.h
+++ b/Solver/functionSpace.h
@@ -12,7 +12,6 @@
#include "simpleFunction.h"
//class SVoid{};
-
template<class T> struct TensorialTraits
{
typedef T ValType;
@@ -141,7 +140,7 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
gradsuvw[i][2]));
}
- virtual int getNumKeys(MElement *ele)
+ virtual int getNumKeys(MElement *ele)
{
if (ele->getParent()) ele = ele->getParent();
return ele->getNumVertices();
diff --git a/Solver/linearSystemCSR.cpp b/Solver/linearSystemCSR.cpp
index a86ef09e1c..8fe68f52e0 100644
--- a/Solver/linearSystemCSR.cpp
+++ b/Solver/linearSystemCSR.cpp
@@ -130,7 +130,7 @@ const unsigned int M_sort2 = 7;
static void free_ivector(int *v, long nl, long nh)
{
- // free an int vector allocated with ivector()
+ // free an int vector allocated with ivector()
free((char*)(v + nl - 1));
}
diff --git a/Solver/linearSystemGMM.h b/Solver/linearSystemGMM.h
index df0df546f8..6198561a6a 100644
--- a/Solver/linearSystemGMM.h
+++ b/Solver/linearSystemGMM.h
@@ -46,6 +46,7 @@ class linearSystemGmm : public linearSystem<scalar> {
}
_a = 0;
}
+
virtual void addToMatrix(int row, int col, const scalar &val)
{
if(val != 0.0) (*_a)(row, col) += val;
@@ -54,10 +55,12 @@ class linearSystemGmm : public linearSystem<scalar> {
{
val = (*_a)(row, col);
}
+
virtual void addToRightHandSide(int row, const scalar &val)
{
if(val != 0.0) (*_b)[row] += val;
}
+
virtual void getFromRightHandSide(int row, scalar &val) const
{
val = (*_b)[row];
@@ -70,14 +73,14 @@ class linearSystemGmm : public linearSystem<scalar> {
{
gmm::clear(*_a);
}
- virtual void zeroRightHandSide()
+ virtual void zeroRightHandSide()
{
for(unsigned int i = 0; i < _b->size(); i++) (*_b)[i] = 0.;
}
void setPrec(double p){ _prec = p; }
void setNoisy(int n){ _noisy = n; }
void setGmres(int n){ _gmres = n; }
- virtual int systemSolve()
+ virtual int systemSolve()
{
//gmm::ilutp_precond<gmm::row_matrix<gmm::wsvector<scalar> > > P(*_a, 25, 0.);
gmm::ildltt_precond<gmm::row_matrix<gmm::wsvector<scalar> > > P(*_a, 30, 1.e-10);
@@ -112,7 +115,7 @@ class linearSystemGmm : public linearSystem<scalar> {
virtual void clear(){}
void setNoisy(int n){}
void setGmres(int n){}
-};
+};
#endif
diff --git a/Solver/terms.h b/Solver/terms.h
index e7d22cb293..f629943724 100644
--- a/Solver/terms.h
+++ b/Solver/terms.h
@@ -27,12 +27,6 @@ class BilinearTermBase
public :
virtual ~BilinearTermBase() {}
virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m) = 0 ;
- /* has nothing to do here...
- virtual void getInter(MInterfaceElement *iele,int npts,IntPt *GP,fullMatrix<double> &m){} ;
- virtual void getInterForce(MInterfaceElement *iele,int npts,IntPt *GP,const fullMatrix<double> &disp,fullMatrix<double> &m){} ;
- virtual void getForce(MElement *ele,int npts,IntPt *GP,const fullMatrix<double> &disp, fullMatrix<double> &m){}
- virtual void getInterOnBoundary(MInterfaceElement *iele,int npts,IntPt *GP,fullMatrix<double> &m){}
- virtual void getInterForceOnBoundary(MInterfaceElement *iele,int npts,IntPt *GP,const fullMatrix<double> &disp,fullMatrix<double> &m){} ;*/
};
template<class T1,class T2> class BilinearTerm : public BilinearTermBase
diff --git a/contrib/cm3/AssembleInterface.h b/contrib/cm3/AssembleInterface.h
deleted file mode 100644
index be1e55d9e8..0000000000
--- a/contrib/cm3/AssembleInterface.h
+++ /dev/null
@@ -1,54 +0,0 @@
-//
-// C++ Interface: terms
-//
-// Description: Function to Assemble the interface term
-// TODO Template this function with Solver/solverAlgorithms
-//
-// Author: <Gauthier BECKER>, (C) 2010
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-
-
-#ifndef ASSEMBLEINTERFACE_H_
-#define ASSEMBLEINTERFACE_H_
-#include "dofManager.h"
-#include "C0DgPlateTerms.h"
-#include "quadratureRules.h"
-#include "MVertex.h"
-#include "MInterfaceElement.h"
-
-void AssembleInterface(BilinearTermBase &iterm,FunctionSpaceBase &space, groupOfElements* ve, std::vector<MInterfaceElement*> &vie, QuadratureBase &integrator,dofManager<double> &assembler){
- // creation of a matrix this matrix is resize and put = 0 in the function iterm.get()
- fullMatrix<double> localMatrix;
- // vector with dof the dof are append in space.getKeys()
- std::vector<Dof> R;
- // Loop on interface element
- for (int i = 0; i < vie.size(); i++)
- {
- MInterfaceElement *ie = vie[i];
- // get pointers on the elements linked to the interface element i
- MElement **elem = ie->getElem();
- // Clear the dof vector (has to be clear because getKeys is an append method)
- R.clear();
- // Compute the integration point
- IntPt *GP;
- int npts=integrator.getIntPoints(ie,&GP);
- // compute elementatry matrix
- if(ie->getElem(1) == ie->getElem(0)){ // Element is on the boundary
- iterm.getInterOnBoundary(ie,npts,GP,localMatrix);
- space.getKeys(ie->getElem(0),R); // dof
- }
- else{
- iterm.getInter(ie,npts,GP,localMatrix);
- // Dof of minus and plus elements
- space.getKeys(ie->getElem(0),R);
- space.getKeys(ie->getElem(1),R);
- }
- assembler.assemble(R, localMatrix);
- }
-
-}
-
-#endif // ASSEMBLEINTERFACE_H_
diff --git a/contrib/cm3/C0DgPlateTerms.h b/contrib/cm3/C0DgPlateTerms.h
deleted file mode 100644
index 5048a2b857..0000000000
--- a/contrib/cm3/C0DgPlateTerms.h
+++ /dev/null
@@ -1,650 +0,0 @@
-//
-// C++ Interface: terms
-//
-// Description: Elementary matrix terms for C0 Dg Plate
-//
-//
-// Author: <Gauthier BECKER>, (C) 2010
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-
-
-
-#ifndef C0DGPLATETERMS_H
-#define C0DGPLATETERMS_H
-#include "terms.h"
-#include "SVector3.h"
-#include <string>
-#include "LocalBasis.h"
-#include "LinearElasticShellHookeTensor.h"
-#include "DgC0PlateElementaryTerms.h"
-
-
-class IsotropicElasticTermC0Plate : public BilinearTerm<SVector3,SVector3>
-{
- protected :
- double E,nu,h,Cm,Cn;
- bool sym;
-
- public :
- IsotropicElasticTermC0Plate(FunctionSpace<SVector3>& space1_,FunctionSpace<SVector3>& space2_,double E_,double nu_,double h_) : BilinearTerm<SVector3,SVector3>(space1_,space2_),E(E_),nu(nu_),h(h_)
- {
- sym=(&space1_==&space2_);
- Cm = ( E_ * h_ * h_ * h_ ) / ( 12 * (1 - nu_) * (1 + nu_) );
- Cn = E_*h_/((1-nu_)*(1+nu_));
- }
- IsotropicElasticTermC0Plate(FunctionSpace<SVector3>& space1_,double E_,double nu_,double h_) : BilinearTerm<SVector3,SVector3>(space1_,space1_),E(E_),nu(nu_),h(h_)
- {
- sym=true;
- Cm = ( E_ * h_ * h_ * h_ ) / ( 12 * (1 - nu_) * (1 + nu_) );
- Cn = E_*h_/((1-nu_)*(1+nu_));
- }
- virtual ~IsotropicElasticTermC0Plate() {}
- virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m)
- {
- if (sym)
- {
- // Initialization of some data
- const int nbdof = BilinearTerm<SVector3,SVector3>::space1.getNumKeys(ele);
- const int nbFF = nbdof/3;
- double Cmt= 0., Cnt=0.;
- LinearElasticShellHookeTensor HOOKe; LinearElasticShellHookeTensor *H=&HOOKe; // make a pointer in an other way
- m.resize(nbdof, nbdof);
- m.setAll(0.);
- std::vector<TensorialTraits<SVector3>::HessType> Hess;
- std::vector<TensorialTraits<SVector3>::GradType> Grad;
- double Bn[256][3][2][2]; // max order 256 or dynamical allocation ?? better than std::vector and fullMatrix ?? create a class with this variable
- LocalBasis LB;
- LocalBasis *lb=&LB; // two last line in 1 operation ??
- double me_bending[3][3];
- double me_membrane[3][3]; // better than fullMatrix<double> used for me_bending ??
- // sum on Gauss' points
- for (int i = 0; i < npts; i++)
- {
- // Coordonate of Gauss' point i
- const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2];
- // Weight of Gauss' point i and Jacobian's value at this point
- const double weight = GP[i].weight;
-
- // Compute of Hessian of SF. Each component of the vector in link to a shape function.
- // It give the value of second derivative of SF in the isoparametric configuration
- BilinearTerm<SVector3,SVector3>::space1.gradfuvw(ele,u, v, w, Grad); // a optimiser the jacobian cannot be given in argument...
- BilinearTerm<SVector3,SVector3>::space1.hessfuvw(ele,u, v, w, Hess); // a optimiser
-
- lb->set(ele,Grad); // This function can become a method of MElement Plante lors du calcul de l'énergie à cause de std::vector<TensorialTraits<SVector3>::GradType> qui retourne un vecteur vide prob de template??
-
- // multiplication of constant by detJ and weight
- Cmt = lb->getJacobian() * weight *Cm;
- Cnt = lb->getJacobian() * weight *Cn;
- // compute of Hooke tensor
- H->set(lb,nu);
-
- // compute Bn value
- compute_Bn(lb,Grad,nbFF,Bn);
-
- // loop on SF to construct the elementary (bulk) stiffness matrix at the Gauss' point i
- for(int j=0; j<nbFF;j++)
- for(int k=0;k<nbFF;k++){
- BulkC0PlateDGStiffnessBendingTerms(Hess[j],Hess[k],H,lb,me_bending);
- BulkC0PlateDGStiffnessMembraneTerms(Bn[j],Bn[k],H,me_membrane);
- for(int jj=0;jj<3;jj++)
- for(int kk=0;kk<3;kk++)
- m(j+(jj*nbFF),k+(kk*nbFF)) += (Cmt*me_bending[jj][kk]+Cnt*me_membrane[jj][kk]);
- }
- // clear the hessian and Grad because the components append in hessfuvw and gradfuvw
- Hess.clear(); Grad.clear();
- }
-// m.print("bulk");
-/* // By numerical perturbation (Verification OK)
- double eps=1.e-8;
- fullMatrix<double> fp(nbdof,1);
- fullMatrix<double> fm(nbdof,1);
- fp.setAll(0.);
- fm.setAll(0.);
- fullMatrix<double> m2(nbdof,nbdof); m2.setAll(0.);
- fullMatrix<double> disp(nbdof,1);
- //GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
- //IntPt *GPbulk;
- //int npts=integrator.getIntPoints(ele,&GP);
- for(int j=0;j<nbdof;j++){
- disp.setAll(0.); fm.setAll(0.);fp.setAll(0.);
- disp(j,0)=eps;
- this->getForce(ele,npts,GP,disp,fp);
- disp(j,0)=-eps;
- this->getForce(ele,npts,GP,disp,fm);
- for(int k=0;k<nbdof;k++)
- m2(k,j)=(fp(k,0)-fm(k,0))/(2.*eps);
- }
- m2.print("Bulk pert");*/
- }
- else
- {
- printf("not implemented\n");
- }
- }
-
-
- virtual void getForce(MElement *ele,int npts,IntPt *GP,const fullMatrix<double> &disp, fullMatrix<double> &m)
- {
- if (sym)
- {
- // Initialization of some data
- const int nbdof = BilinearTerm<SVector3,SVector3>::space1.getNumKeys(ele);
- const int nbFF = nbdof/3;
- double Cmt= 0.,Cnt=0.;
- LinearElasticShellHookeTensor HOOKe; LinearElasticShellHookeTensor *H=&HOOKe; // make a pointer in an other way
- m.resize(nbdof,1);
- m.setAll(0.);
- std::vector<TensorialTraits<SVector3>::HessType> Hess;
- std::vector<TensorialTraits<SVector3>::GradType> Grad;
- double Bn[256][3][2][2]; // max order 256 or dynamical allocation ?? better than std::vector and fullMatrix ?? create a class with this variable
- LocalBasis LB;
- LocalBasis *lb=&LB; // two last line in 1 operation ??
- double me_bending[3];
- double me_membrane[3];
- // sum on Gauss' points
- for (int i = 0; i < npts; i++)
- {
- // Coordonate of Gauss' point i
- const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2];
- // Weight of Gauss' point i and Jacobian's value at this point
- const double weight = GP[i].weight;
-
- // Compute of Hessian of SF. Each component of the vector in link to a shape function.
- // It give the value of second derivative of SF in the isoparametric configuration
- BilinearTerm<SVector3,SVector3>::space1.gradfuvw(ele,u, v, w, Grad); // a optimiser the jacobian cannot be given in argument...
- BilinearTerm<SVector3,SVector3>::space1.hessfuvw(ele,u, v, w, Hess); // a optimiser
-
- lb->set(ele,Grad); // This function can become a method of MElement Plante lors du calcul de l'énergie à cause de std::vector<TensorialTraits<SVector3>::GradType> qui retourne un vecteur vide prob de template??
-
- // multiplication of constant by detJ and weight
- Cmt = lb->getJacobian() * weight * Cm;
- Cnt = lb->getJacobian() * weight * Cn;
-
- // compute of Hooke tensor
- H->set(lb,nu);
-
- // compute Bn value
- compute_Bn(lb,Grad,nbFF,Bn);
-
- // loop on SF to construct the elementary (bulk) stiffness matrix at the Gauss' point i
- for(int j=0; j<nbFF;j++){
- BulkC0PlateDGForceBendingTerms(Hess[j],Hess,H,lb,disp,me_bending);
- BulkC0PlateDGForceMembraneTerms(Bn[j],Bn,H,nbFF,disp,me_membrane);
- for(int k=0;k<3;k++)
- m(j+k*nbFF,0) += (Cmt*me_bending[k]+ Cnt*me_membrane[k]);
- }
- // clear the hessian and Grad because the components append in hessfuvw and gradfuvw
- Hess.clear(); Grad.clear();
- }
- }
- else
- {
- printf("not implemented\n");
- }
- }
-}; // IsotropicElasticTermC0Plate
-
-
-class IsotropicElasticInterfaceTermC0Plate : public BilinearTerm<SVector3,SVector3>
-{
- protected :
- double E,nu,beta1,h;
- double Cm;
- bool sym;
-
- public :
- IsotropicElasticInterfaceTermC0Plate(FunctionSpace<SVector3>& space1_,FunctionSpace<SVector3>& space2_,double E_,double nu_,double beta1_,double h_) : BilinearTerm<SVector3,SVector3>(space1_,space2_),E(E_),nu(nu_),beta1(beta1_),h(h_)
- {
- Cm = ( E_ * h_ * h_ * h_ ) / ( 12 * (1 - nu_) * (1 + nu_) );
- sym=(&space1_==&space2_);
- }
- IsotropicElasticInterfaceTermC0Plate(FunctionSpace<SVector3>& space1_,double E_,double nu_,double beta1_,double h_) : BilinearTerm<SVector3,SVector3>(space1_,space1_),E(E_),nu(nu_),beta1(beta1_),h(h_)
- {
- Cm = ( E_ * h_ * h_ * h_ ) / ( 12 * (1 - nu_) * (1 + nu_) );
- sym=true;
- }
- virtual ~IsotropicElasticInterfaceTermC0Plate() {}
- // Must be declare has this function is virtual pure in class BilinearTermBase
- virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m) {}
-
- virtual void getInter(MInterfaceElement *iele,int npts,IntPt *GP,fullMatrix<double> &m)
- {
- if (sym)
- {
- // data initialization
- // Retrieve of the element link to the interface element velem[0] == minus elem ; velem == plus elem
- MElement ** velem = iele->getElem();
- const int nbdof_m = BilinearTerm<SVector3,SVector3>::space1.getNumKeys(velem[0]);
- const int nbdof_p = BilinearTerm<SVector3,SVector3>::space1.getNumKeys(velem[1]);
- const int nbFF_m = nbdof_m/3;
- const int nbFF_p = nbdof_p/3;
- // Initialization
- m.resize(nbdof_m+nbdof_p, nbdof_m+nbdof_p);
- m.setAll(0.);
- double uem,uep,vem,vep;
- std::vector<TensorialTraits<SVector3>::GradType> Grads_m;
- std::vector<TensorialTraits<SVector3>::GradType> Grads_p;
- std::vector<TensorialTraits<SVector3>::GradType> Grads;
- std::vector<TensorialTraits<SVector3>::HessType> Hess_m;
- std::vector<TensorialTraits<SVector3>::HessType> Hess_p;
- LocalBasis LBS,LBP,LBM;
- LocalBasis *lbs=&LBS; LocalBasis *lb_p=&LBP; LocalBasis *lb_m=&LBM;
- double Bhat_p[256][3][2][2], Bhat_m[256][3][2][2];
- LinearElasticShellHookeTensor HOOKEhat_p;
- LinearElasticShellHookeTensor HOOKehat_m;
- LinearElasticShellHookeTensor HOOKehatmean;
- LinearElasticShellHookeTensor *Hhat_p=&HOOKEhat_p;
- LinearElasticShellHookeTensor *Hhat_m=&HOOKehat_m;
- LinearElasticShellHookeTensor *Hhatmean=&HOOKehatmean;
- double Deltat_m[256][3][3], Deltat_p[256][3][3];
- double Cmt;
- double me_cons[3][3];
- double me_comp[3][3];
- double me_stab[3][3];
-
- // Characteristic size of interface element
- double h_s = iele->getCharacteristicSize();
-
- const double Bhs = beta1/h_s;
- // sum on Gauss' points
- for (int i = 0; i < npts; i++)
- {
- // Coordonate of Gauss' point i
- const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2];
- // Weight of Gauss' point i
- const double weight = GP[i].weight;
- //printf("Abscisse of gauss point %f\n",u);
- //Compute the position (u,v) in the element of the Gauss point (to know where evaluate the shape functions)
- iele->getuvOnElem(u,uem,vem,uep,vep);
- // Compute of gradient and hessian of shape functions on interface element and on elements
- // ATTENTION after multiplication by multipliers (functionspace 276) components are in the "second line of tensor change this ??
- BilinearTerm<SVector3,SVector3>::space1.gradfuvw(iele,u, v, w, Grads); // grad of shape fonction on interface element
- BilinearTerm<SVector3,SVector3>::space1.gradfuvw(velem[0],uem, vem, w, Grads_m); // w = 0
- BilinearTerm<SVector3,SVector3>::space1.gradfuvw(velem[1],uep, vep, w, Grads_p);
- BilinearTerm<SVector3,SVector3>::space1.hessfuvw(velem[0],uem, vem, w, Hess_m);
- BilinearTerm<SVector3,SVector3>::space1.hessfuvw(velem[1],uep, vep, w, Hess_p);
-
- // basis of elements and interface element
- lb_m->set(velem[0],Grads_m); // This function can become a method of MElement Plante lors du calcul de l'énergie à cause de std::vector<TensorialTraits<SVector3>::GradType> prob de template??
- lb_p->set(velem[1],Grads_p); // This function can become a method of MElement Plante lors du calcul de l'énergie à cause de std::vector<TensorialTraits<SVector3>::GradType> prob de template??
- lbs->set(iele,Grads,lb_p->gett0(),lb_m->gett0()); // This function can become a method of MElement Plante lors du calcul de l'énergie à cause de std::vector<TensorialTraits<SVector3>::GradType> prob de template??
-
- // PushForwardTensor
- lb_m->set_pushForward(lbs);
- lb_p->set_pushForward(lbs);
-
- // Compute of Bhat vector (1 component for now because 1 dof (z) ) --> it's a vector of length == nbFF
- Compute_Bhat(lb_p,Hess_p,nbFF_p,Bhat_p);
- Compute_Bhat(lb_m,Hess_m,nbFF_m,Bhat_m);
-
- // Compute of Hooke hat tensor on minus and plus element
- Cmt = weight * lbs->getJacobian() * Cm; // Eh^3/(12(1-nu^2)) * weight gauss * jacobian
- Hhat_p->hat(lb_p,Cmt,nu);
- Hhat_m->hat(lb_m,Cmt,nu); //Hhat_m->hat(lb_m,H_m);
- Hhatmean->mean(Hhat_p,Hhat_m); // mean value of tensor by component used for stability term
-
- // Compute of Deltat_tilde
- compute_Deltat_tilde(lb_p,Grads_p,nbFF_p,Deltat_p);
- compute_Deltat_tilde(lb_m,Grads_m,nbFF_m,Deltat_m);
-
- for(int j=0;j<nbFF_m;j++){
- for(int k=0;k<nbFF_m;k++){
- consAndCompC0PlateStiffnessTerms(Hhat_m,Bhat_m[j],Deltat_m[k],lbs,me_cons);
- consAndCompC0PlateStiffnessTerms(Hhat_m,Bhat_m[k],Deltat_m[j],lbs,me_comp);
- stabilityC0PlateStiffnessTerms(Hhatmean,Deltat_m[j], Deltat_m[k],lbs,me_stab);
- for(int jj=0;jj<3;jj++)
- for(int kk=0;kk<3;kk++)
- m(j+jj*nbFF_m,k+kk*nbFF_m) += (- me_cons[jj][kk] - me_comp[jj][kk] + Bhs * me_stab[jj][kk] );
- }
- for(int k=0;k<nbFF_p;k++){
- consAndCompC0PlateStiffnessTerms(Hhat_m,Bhat_m[j],Deltat_p[k],lbs,me_cons);
- consAndCompC0PlateStiffnessTerms(Hhat_p,Bhat_p[k],Deltat_m[j],lbs,me_comp);
- stabilityC0PlateStiffnessTerms(Hhatmean,Deltat_m[j], Deltat_p[k],lbs,me_stab);
- for(int jj=0;jj<3;jj++)
- for(int kk=0;kk<3;kk++)
- m(j+jj*nbFF_m,k+nbdof_m+kk*nbFF_p) += ( me_cons[jj][kk] - me_comp[jj][kk]- Bhs * me_stab[jj][kk] );
- }
- }
- for(int j=0;j<nbFF_p;j++){
- for(int k=0;k<nbFF_m;k++){
- consAndCompC0PlateStiffnessTerms(Hhat_p,Bhat_p[j],Deltat_m[k],lbs,me_cons);
- consAndCompC0PlateStiffnessTerms(Hhat_m,Bhat_m[k],Deltat_p[j],lbs,me_comp);
- stabilityC0PlateStiffnessTerms(Hhatmean,Deltat_p[j], Deltat_m[k],lbs,me_stab);
- for(int jj=0;jj<3;jj++)
- for(int kk=0;kk<3;kk++)
- m(j+(jj*nbFF_p)+nbdof_m,k+kk*nbFF_m) += (- me_cons[jj][kk] + me_comp[jj][kk] - Bhs * me_stab[jj][kk]);
- }
- for(int k=0;k<nbFF_p;k++){
- consAndCompC0PlateStiffnessTerms(Hhat_p,Bhat_p[j],Deltat_p[k],lbs,me_cons);
- consAndCompC0PlateStiffnessTerms(Hhat_p,Bhat_p[k],Deltat_p[j],lbs,me_comp);
- stabilityC0PlateStiffnessTerms(Hhatmean,Deltat_p[j], Deltat_p[k],lbs,me_stab);
- for(int jj=0;jj<3;jj++)
- for(int kk=0;kk<3;kk++)
- m(j+(jj*nbFF_p)+nbdof_m,k+(kk*nbFF_p)+nbdof_m) += ( me_cons[jj][kk] + me_comp[jj][kk] + Bhs * me_stab[jj][kk]);
- }
- }
- // Because component are push_back in Grads in gradfuvw idem for hess
- Grads_m.clear(); Grads_p.clear(); Hess_m.clear(); Hess_p.clear(); Grads.clear();
- }
-// m.print("Interface");
-/* // Compute the matrix by numerical perturbation Verif OK
- double eps=1.e-6;
- fullMatrix<double> fp(nbdof_m+nbdof_p,1);
- fullMatrix<double> fm(nbdof_m+nbdof_p,1);
- fp.setAll(0.);
- fm.setAll(0.);
- fullMatrix<double> m2(nbdof_m+nbdof_p,nbdof_m+nbdof_p); m2.setAll(0.);
- fullMatrix<double> disp(nbdof_m+nbdof_p,1);
- for(int j=0;j<nbdof_m+nbdof_p;j++){
- disp.setAll(0.); fm.setAll(0.);fp.setAll(0.);
- disp(j,0)=eps;
- this->getInterForce(iele,npts,GP,disp,fp);
- disp(j,0)=-eps;
- this->getInterForce(iele,npts,GP,disp,fm);
- for(int k=0;k<nbdof_m+nbdof_p;k++)
- m2(k,j)=(fp(k,0)-fm(k,0))/(2.*eps);
- }
- m2.print("Matrix by perturbation");*/
- }
- else
- {
- printf("not implemented\n");
- }
- }
-
- virtual void getInterForce(MInterfaceElement *iele,int npts,IntPt *GP,const fullMatrix<double> &disp, fullMatrix<double> &m){
- if (sym)
- {
- // data initialization
- // Retrieve of the element link to the interface element velem[0] == minus elem ; velem == plus elem
- MElement ** velem = iele->getElem();
- const int nbdof_m = BilinearTerm<SVector3,SVector3>::space1.getNumKeys(velem[0]);
- const int nbdof_p = BilinearTerm<SVector3,SVector3>::space1.getNumKeys(velem[1]);
- const int nbFF_m=nbdof_m/3;
- const int nbFF_p=nbdof_p/3;
- // Initialization
- m.resize(nbdof_m+nbdof_p, 1);
- m.setAll(0.);
- double uem,uep,vem,vep;
- std::vector<TensorialTraits<SVector3>::GradType> Grads_m;
- std::vector<TensorialTraits<SVector3>::GradType> Grads_p;
- std::vector<TensorialTraits<SVector3>::GradType> Grads;
- std::vector<TensorialTraits<SVector3>::HessType> Hess_m;
- std::vector<TensorialTraits<SVector3>::HessType> Hess_p;
- LocalBasis LBS,LBP,LBM;
- LocalBasis *lbs=&LBS; LocalBasis *lb_p=&LBP; LocalBasis *lb_m=&LBM;
- double Bhat_p[256][3][2][2],Bhat_m[256][3][2][2];
- LinearElasticShellHookeTensor HOOKEhat_p; LinearElasticShellHookeTensor *Hhat_p=&HOOKEhat_p;
- LinearElasticShellHookeTensor HOOKehat_m; LinearElasticShellHookeTensor *Hhat_m=&HOOKehat_m;
- LinearElasticShellHookeTensor HOOKehatmean; LinearElasticShellHookeTensor *Hhatmean=&HOOKehatmean;
- double Deltat_m[256][3][3], Deltat_p[256][3][3];
- double Cmt;
- double me_cons[3];
- double me_comp[3];
- double me_stab[3];
-
- // Characteristic size of interface element
- double h_s = iele->getCharacteristicSize();
-
- const double Bhs = beta1/h_s;
- // sum on Gauss' points
- for (int i = 0; i < npts; i++)
- {
- // Coordonate of Gauss' point i
- const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2];
- // Weight of Gauss' point i
- const double weight = GP[i].weight;
- //printf("Abscisse of gauss point %f\n",u);
- //Compute the position (u,v) in the element of the Gauss point (to know where evaluate the shape functions)
- iele->getuvOnElem(u,uem,vem,uep,vep);
- // Compute of gradient and hessian of shape functions on interface element and on elements
- // ATTENTION after multiplication by multipliers (functionspace 276) components are in the "second line of tensor change this ??
- BilinearTerm<SVector3,SVector3>::space1.gradfuvw(iele,u, v, w, Grads); // grad of shape fonction on interface element
- BilinearTerm<SVector3,SVector3>::space1.gradfuvw(velem[0],uem, vem, w, Grads_m); // w = 0
- BilinearTerm<SVector3,SVector3>::space1.gradfuvw(velem[1],uep, vep, w, Grads_p);
- BilinearTerm<SVector3,SVector3>::space1.hessfuvw(velem[0],uem, vem, w, Hess_m);
- BilinearTerm<SVector3,SVector3>::space1.hessfuvw(velem[1],uep, vep, w, Hess_p);
-
- // basis of elements and interface element
- lb_m->set(velem[0],Grads_m); // This function can become a method of MElement Plante lors du calcul de l'énergie à cause de std::vector<TensorialTraits<SVector3>::GradType> prob de template??
- lb_p->set(velem[1],Grads_p); // This function can become a method of MElement Plante lors du calcul de l'énergie à cause de std::vector<TensorialTraits<SVector3>::GradType> prob de template??
- lbs->set(iele,Grads,lb_p->gett0(),lb_m->gett0()); // This function can become a method of MElement Plante lors du calcul de l'énergie à cause de std::vector<TensorialTraits<SVector3>::GradType> prob de template??
-
- // PushForwardTensor
- lb_m->set_pushForward(lbs);
- lb_p->set_pushForward(lbs);
-
- // Compute of Bhat vector (1 component for now because 1 dof (z) ) --> it's a vector of length == nbFF
- Compute_Bhat(lb_p,Hess_p,nbFF_p,Bhat_p);
- Compute_Bhat(lb_m,Hess_m,nbFF_m,Bhat_m);
- // Compute of Hooke hat tensor on minus and plus element
- Cmt = weight * lbs->getJacobian() * Cm; // Eh^3/(12(1-nu^2)) * weight gauss * jacobian
- Hhat_p->hat(lb_p,Cmt,nu);
- Hhat_m->hat(lb_m,Cmt,nu);
- Hhatmean->mean(Hhat_p,Hhat_m); // mean value of tensor by component used for stability term
-
- // Compute of Deltat_tilde
- compute_Deltat_tilde(lb_p,Grads_p,nbFF_p,Deltat_p);
- compute_Deltat_tilde(lb_m,Grads_m,nbFF_m,Deltat_m);
-
- // loop on dof ATTENTION SAME NUMBER OF DOF for the two elements TODO take into account a different dof numbers between the two elements There ok because sym ??
- for(int j=0;j<nbFF_m;j++){
- consC0PlateForceTerms(Hhat_m,nbFF_m,nbFF_p,Bhat_m[j],Deltat_m,Deltat_p,lbs,disp,me_cons);
- compC0PlateForceTerms(nbFF_m,nbFF_p, Hhat_m,Hhat_p,Bhat_m,Bhat_p,Deltat_m[j],lbs,disp,me_comp);
- stabilityC0PlateForceTerms(nbFF_m,nbFF_p,Hhatmean,Deltat_m[j],Deltat_m,Deltat_p,lbs,disp,me_stab);
- for(int jj=0;jj<3;jj++)
- m(j+jj*nbFF_m,0) += (me_cons[jj] - me_comp[jj]- Bhs * me_stab[jj] );
- }
- for(int j=0;j<nbFF_p;j++){
- consC0PlateForceTerms(Hhat_p,nbFF_m,nbFF_p,Bhat_p[j],Deltat_m,Deltat_p,lbs,disp,me_cons);
- compC0PlateForceTerms(nbFF_m,nbFF_p,Hhat_m,Hhat_p,Bhat_m,Bhat_p,Deltat_p[j],lbs,disp,me_comp);
- stabilityC0PlateForceTerms(nbFF_m,nbFF_p,Hhatmean,Deltat_p[j],Deltat_m,Deltat_p,lbs,disp,me_stab);
- for(int jj=0;jj<3;jj++)
- m(j+(jj*nbFF_p)+nbdof_m,0) += ( me_cons[jj] + me_comp[jj] + Bhs * me_stab[jj]);
- }
- // Because component are push_back in Grads in gradfuvw idem for hess
- Grads_m.clear(); Grads_p.clear(); Hess_m.clear(); Hess_p.clear(); Grads.clear();
- }
- }
- else
- {
- printf("not implemented\n");
- }
-
- }
-
-
- virtual void getInterOnBoundary(MInterfaceElement *iele,int npts,IntPt *GP,fullMatrix<double> &m)
- {
- if (sym)
- {
- // data initialization
- // Retrieve of the element link to the interface element velem[0] == minus elem ; velem == plus elem
- MElement ** velem = iele->getElem();
- const int nbdof_m = BilinearTerm<SVector3,SVector3>::space1.getNumKeys(velem[0]);
- const int nbFF_m = nbdof_m/3;
- // Initialization
- m.resize(nbdof_m, nbdof_m);
- m.setAll(0.);
- double uem,uep,vem,vep;
- std::vector<TensorialTraits<SVector3>::GradType> Grads_m;
- std::vector<TensorialTraits<SVector3>::GradType> Grads_p;
- std::vector<TensorialTraits<SVector3>::GradType> Grads;
- std::vector<TensorialTraits<SVector3>::HessType> Hess_m;
- std::vector<TensorialTraits<SVector3>::HessType> Hess_p;
- LocalBasis LBS,LBM;
- LocalBasis *lbs=&LBS; LocalBasis *lb_m=&LBM;
- //std::vector<std::vector<fullMatrix<double> > > Bhat_m;
- double Bhat_m[256][3][2][2];
- LinearElasticShellHookeTensor HOOKehat_m;
- LinearElasticShellHookeTensor *Hhat_m=&HOOKehat_m;
- double Deltat_m[256][3][3],Deltat_p[256][3][3];
- double Cmt;
- double me_cons[3][3];
- double me_comp[3][3];
- double me_stab[3][3];
-
- // Characteristic size of interface element
- double h_s = iele->getCharacteristicSize();
- const double Bhs = beta1/h_s;
- // sum on Gauss' points
- for (int i = 0; i < npts; i++)
- {
- // Coordonate of Gauss' point i
- const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2];
- // Weight of Gauss' point i
- const double weight = GP[i].weight;
-
- //printf("Abscisse of gauss point %f\n",u);
- //Compute the position (u,v) in the element of the Gauss point (to know where evaluate the shape functions)
- iele->getuvOnElem(u,uem,vem,uep,vep);
- //printf("Position (u,v) of minus element (%f,%f)\n",uem,vem);
- // Compute of gradient and hessian of shape functions on interface element and on elements
- // ATTENTION after multiplication by multipliers (functionspace 276) components are in the "second line of tensor change this ??
- BilinearTerm<SVector3,SVector3>::space1.gradfuvw(iele,u, v, w, Grads); // grad of shape fonction on interface element
- BilinearTerm<SVector3,SVector3>::space1.gradfuvw(velem[0],uem, vem, w, Grads_m); // w = 0
- BilinearTerm<SVector3,SVector3>::space1.hessfuvw(velem[0],uem, vem, w, Hess_m);
-
- // basis of elements and interface element
- lb_m->set(velem[0],Grads_m); // This function can become a method of MElement Plante lors du calcul de l'énergie à cause de std::vector<TensorialTraits<SVector3>::GradType> prob de template??
- lbs->set(iele,Grads,lb_m->gett0()); // This function can become a method of MElement Plante lors du calcul de l'énergie à cause de std::vector<TensorialTraits<SVector3>::GradType> prob de template??
-
- // PushForwardTensor
- lb_m->set_pushForward(lbs);
-
- /*printf("Local basis interface\n");
- lbs->print();
- printf("Local basis minus\n");
- lb_m->print();*/
-
- // Compute of Bhat vector (1 component for now because 1 dof (z) ) --> it's a vector of length == nbFF
- Compute_Bhat(lb_m,Hess_m,nbFF_m,Bhat_m);
- // Compute of Hooke hat tensor on minus and plus element
- Cmt = weight * lbs->getJacobian() * Cm; // Eh^3/(12(1-nu^2)) * weight gauss * jacobian
- Hhat_m->hat(lb_m,Cmt,nu);
-
- // Compute of Deltat_tilde
- compute_Deltat_tildeBound(lb_m,Grads_m,nbFF_m,Deltat_m,lbs);
-
- // loop on dof ATTENTION SAME NUMBER OF DOF for the two elements TODO take into account a different dof numbers between the two elements There ok because sym ??
- for(int j=0;j<nbFF_m;j++){
- for(int k=0;k<nbFF_m;k++){
- consAndCompC0PlateStiffnessTerms(Hhat_m,Bhat_m[j],Deltat_m[k],lbs,me_cons);
- consAndCompC0PlateStiffnessTerms(Hhat_m,Bhat_m[k],Deltat_m[j],lbs,me_comp);
- stabilityC0PlateStiffnessTerms(Hhat_m,Deltat_m[j], Deltat_m[k],lbs,me_stab);
- for(int jj=0;jj<3;jj++)
- for(int kk=0;kk<3;kk++)
- m(j+jj*nbFF_m,k+kk*nbFF_m) += -(me_cons[jj][kk] + me_comp[jj][kk] - Bhs * me_stab[jj][kk]);
- }
- }
-
- // Because component are push_back in Grads in gradfuvw idem for hess
- Grads_m.clear(); Hess_m.clear(); Grads.clear();
- }
-// m.print("InterfaceBound");
-/* // Compute the matrix by numerical perturbation Verif OK
- double eps=1.e-8;
- fullMatrix<double> fm(nbdof_m,1), fp(nbdof_m,1);
- fullMatrix<double> m2(nbdof_m,nbdof_m); m2.setAll(0.);
- fullMatrix<double> disp(nbdof_m,1);
- for(int j=0;j<nbdof_m;j++){
- disp.setAll(0.); fm.setAll(0.); fp.setAll(0.);
- disp(j,0)=eps;
- this->getInterForceOnBoundary(iele,npts,GP,disp,fp);
- disp(j,0)=-eps;
- this->getInterForceOnBoundary(iele,npts,GP,disp,fm);
- for(int k=0;k<nbdof_m;k++)
- m2(k,j)=(fp(k,0)-fm(k,0))/(2.*eps);
- }
- m2.print("Matrix by perturbation");*/
- }
- else
- {
- printf("not implemented\n");
- }
- }
-
- virtual void getInterForceOnBoundary(MInterfaceElement *iele,int npts,IntPt *GP,const fullMatrix<double> &disp, fullMatrix<double> &m){
- if (sym)
- {
- // data initialization
- // Retrieve of the element link to the interface element velem[0] == minus elem ; velem == plus elem
- MElement ** velem = iele->getElem();
- const int nbdof_m = BilinearTerm<SVector3,SVector3>::space1.getNumKeys(velem[0]);
- const int nbFF_m = nbdof_m/3;
- // Initialization
- m.resize(nbdof_m, 1);
- m.setAll(0.);
- double uem,vem,uep,vep;
- std::vector<TensorialTraits<SVector3>::GradType> Grads_m;
- std::vector<TensorialTraits<SVector3>::GradType> Grads;
- std::vector<TensorialTraits<SVector3>::HessType> Hess_m;
- LocalBasis LBS,LBM;
- LocalBasis *lbs=&LBS; LocalBasis *lb_m=&LBM;
- double Bhat_m[256][3][2][2];
- LinearElasticShellHookeTensor HOOKehat_m;
- LinearElasticShellHookeTensor HOOKe_m;
- LinearElasticShellHookeTensor *Hhat_m=&HOOKehat_m;
- double Deltat_m[256][3][3];
- double Cmt;
- double me_cons[3];
- double me_comp[3];
- double me_stab[3];
- // Characteristic size of interface element
- double h_s = iele->getCharacteristicSize();
-
- const double Bhs = beta1/h_s;
- // sum on Gauss' points
- for (int i = 0; i < npts; i++)
- {
- // Coordonate of Gauss' point i
- const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2];
- // Weight of Gauss' point i
- const double weight = GP[i].weight;
- //printf("Abscisse of gauss point %f\n",u);
- //Compute the position (u,v) in the element of the Gauss point (to know where evaluate the shape functions)
- iele->getuvOnElem(u,uem,vem,uep,vep);
- // Compute of gradient and hessian of shape functions on interface element and on elements
- // ATTENTION after multiplication by multipliers (functionspace 276) components are in the "second line of tensor change this ??
- BilinearTerm<SVector3,SVector3>::space1.gradfuvw(iele,u, v, w, Grads); // grad of shape fonction on interface element
- BilinearTerm<SVector3,SVector3>::space1.gradfuvw(velem[0],uem, vem, w, Grads_m); // w = 0
- BilinearTerm<SVector3,SVector3>::space1.hessfuvw(velem[0],uem, vem, w, Hess_m);
-
- // basis of elements and interface element
- lb_m->set(velem[0],Grads_m); // This function can become a method of MElement Plante lors du calcul de l'énergie à cause de std::vector<TensorialTraits<SVector3>::GradType> prob de template??
- lbs->set(iele,Grads,lb_m->gett0()); // This function can become a method of MElement Plante lors du calcul de l'énergie à cause de std::vector<TensorialTraits<SVector3>::GradType> prob de template??
-
- // PushForwardTensor
- lb_m->set_pushForward(lbs);
-
- // Compute of Bhat vector (1 component for now because 1 dof (z) ) --> it's a vector of length == nbFF
- Compute_Bhat(lb_m,Hess_m,nbFF_m,Bhat_m);
- // Compute of Hooke hat tensor on minus and plus element
- Cmt = weight * lbs->getJacobian() * Cm; // Eh^3/(12(1-nu^2)) * weight gauss * jacobian
- Hhat_m->hat(lb_m,Cmt,nu);
-
- // Compute of Deltat_tilde
- compute_Deltat_tildeBound(lb_m,Grads_m,nbFF_m,Deltat_m,lbs);
-
- for(int j=0;j<nbFF_m;j++){
- consC0PlateForceTermsBound(Hhat_m,nbFF_m,Bhat_m[j],Deltat_m,lbs,disp,me_cons);
- compC0PlateForceTermsBound(nbFF_m,Hhat_m,Bhat_m,Deltat_m[j],lbs,disp,me_comp);
- stabilityC0PlateForceTermsBound(nbFF_m,Hhat_m,Deltat_m[j],Deltat_m,lbs,disp,me_stab);
- for(int jj=0;jj<3;jj++)
- m(j+jj*nbFF_m,0) += (me_cons[jj] - me_comp[jj] - Bhs * me_stab[jj] );
- }
- // Because component are push_back in Grads in gradfuvw idem for hess
- Grads_m.clear(); Hess_m.clear(); Grads.clear();
- }
- }
- else
- {
- printf("not implemented\n");
- }
-
- }
-}; // class IsotropicElasticInterfaceTermC0Plate
-#endif // C0DGPLATETERMS_H
diff --git a/contrib/cm3/CMakeLists.txt b/contrib/cm3/CMakeLists.txt
deleted file mode 100644
index 8913c4390a..0000000000
--- a/contrib/cm3/CMakeLists.txt
+++ /dev/null
@@ -1,8 +0,0 @@
-# pour boris:
-## boris cela ne marche que si le fichier est dans SVN ...
-#include(/home/boris/MyGmshProjects/CMakeLists.txt)
-# tests pour eric et christophe:
-list(APPEND EXTERNAL_INCLUDES contrib/cm3)
-add_executable(dgsolver EXCLUDE_FROM_ALL contrib/cm3/mainDG.cpp contrib/cm3/DgC0PlateSolver.cpp contrib/cm3/MInterfaceElement.cpp contrib/cm3/GModelWithInterface.cpp contrib/cm3/C0DgPlateTerms.h ${GMSH_SRC})
-target_link_libraries(dgsolver ${LINK_LIBRARIES})
-
diff --git a/contrib/cm3/DgC0PlateElementaryTerms.h b/contrib/cm3/DgC0PlateElementaryTerms.h
deleted file mode 100644
index b79a215e99..0000000000
--- a/contrib/cm3/DgC0PlateElementaryTerms.h
+++ /dev/null
@@ -1,387 +0,0 @@
-//
-// C++ Interface: terms
-//
-// Description: Functions used to compute a component of elementary matrix term
-//
-//
-// Author: <Gauthier BECKER>, (C) 2010
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-
-// Or put in the file C0DgPlateTerms.h ??
-#ifndef DGC0PLATEELEMENTARYTERMS_H_
-#define DGC0PLATEELEMENTARYTERMS_H_
-
-inline void diaprod(const double a[3], const double b[3], double m[3][3]){
- for(int i=0;i<3;i++)
- for(int j=0;j<3;j++)
- m[i][j]=a[i]*b[j];
-}
-
-static inline void dot(const double a[3], const double b[3], double c[3]){
- for(int i=0;i<3;i++)
- c[i] = a[i]*b[i];
-}
-
-static inline void dot(const double a[3], const SVector3 &b, double c[3]){
- for(int i=0;i<3;i++)
- c[i] = a[i]*b(i);
-}
-
-inline void matTvectprod(const double m[3][3], const SVector3 &v, double v2[3]){
- double temp[3];
- for(int i=0;i<3;i++){
- temp[0] = m[0][i]; temp[1] = m[1][i]; temp[2] = m[2][i];
- v2[i] = dot(temp,v);
- }
-}
-
-inline void matTvectprod(const double m[3][3], const double v[3], double v2[3]){
- double temp[3];
- for(int i=0;i<3;i++){
- temp[0] = m[0][i]; temp[1] = m[1][i]; temp[2] = m[2][i];
- v2[i] = dot(temp,v);
- }
-}
-
-// Compute the component (j,k) of the elementary stiffness matrix
-inline void BulkC0PlateDGStiffnessMembraneTerms(const double Bj[3][2][2],const double Bk[3][2][2], const LinearElasticShellHookeTensor *H, double me[3][3]){
- for(int i=0;i<3;i++)
- for(int j=0;j<3;j++)
- me[i][j]=0.; // make a methode in object me ??
- double mtemp[3][3];
- double v1[3], v2[3];
- for(int a=0;a<2;a++)
- for(int b=0;b<2;b++){
- v2[0] = Bk[0][a][b]; v2[1] = Bk[1][a][b]; v2[2]= Bk[2][a][b]; // make a method
- for(int c=0;c<2;c++)
- for(int d=0;d<2;d++){
- v1[0] = Bj[0][c][d]; v1[1]=Bj[1][c][d]; v1[2]=Bj[2][c][d]; //make a method
- diaprod(v1,v2,mtemp);
- for(int jj=0;jj<3;jj++)
- for(int kk=0;kk<3;kk++)
- me[jj][kk]+=(H->get(a,b,c,d)*mtemp[jj][kk]);
- }
- }
-}
-
-inline void BulkC0PlateDGStiffnessBendingTerms(TensorialTraits<double>::HessType &hessj, TensorialTraits<double>::HessType &hessk, const LinearElasticShellHookeTensor *H, const LocalBasis *lb, double me[3][3]){
- double mtemp[3][3];
- for(int i=0;i<3;i++) for(int j=0;j<3;j++) me[i][j]=0.;
- double B1[3],B2[3];
- for(int alpha=0;alpha<2;alpha++)
- for(int beta=0;beta<2;beta++){
- B1[0]=-hessj(alpha,beta)*lb->gett0(0); B1[1]=-hessj(alpha,beta)*lb->gett0(1); B1[2]=-hessj(alpha,beta)*lb->gett0(2);
- for(int gamma=0;gamma<2;gamma++)
- for(int delta=0;delta<2;delta++){
- B2[0]=-hessk(gamma,delta)*lb->gett0(0); B2[1]=-hessk(gamma,delta)*lb->gett0(1); B2[2]=-hessk(gamma,delta)*lb->gett0(2);
- diaprod(B1,B2,mtemp);
- for(int jj=0;jj<3;jj++)
- for(int kk=0;kk<3;kk++)
- me[jj][kk]+=H->get(alpha,beta,gamma,delta)*mtemp[jj][kk];
- }
- }
-}
-
-inline void consAndCompC0PlateStiffnessTerms(LinearElasticShellHookeTensor *Hhat,const double Bhat[3][2][2],const double dt[3][3], const LocalBasis *lb, double me[3][3]){
- for(int i=0;i<3;i++) for(int j=0;j<3;j++) me[i][j]=0.;
- double v1[3], v2[3];
- double mtemp[3][3];
- double temp=0.;
- for(int alpha=0;alpha<2;alpha++)
- for(int beta=0;beta<2;beta++)
- for(int gamma=0;gamma<2;gamma++)
- for(int delta=0;delta<2;delta++){
- v1[0] = Bhat[0][gamma][delta]; v1[1] = Bhat[1][gamma][delta]; v1[2] = Bhat[2][gamma][delta];
- matTvectprod(dt,lb->getphi0(alpha),v2);
- diaprod(v1,v2,mtemp);
- temp = 0.5*Hhat->get(alpha,beta,gamma,delta)*(-lb->getphi0(1,beta));
- for(int i=0;i<3;i++)
- for(int j=0;j<3;j++)
- me[i][j] += (mtemp[i][j]*temp);
- }
-}
-
-inline void stabilityC0PlateStiffnessTerms(LinearElasticShellHookeTensor *Hhat, const double dta[3][3], const double dtb[3][3], const LocalBasis *lb, double me[3][3]){
- for(int i=0;i<3;i++) for(int j=0;j<3;j++) me[i][j]=0.;
- double mtemp[3][3];
- double v1[3],v2[3];
- double temp=0.;
- for(int alpha=0;alpha<2;alpha++)
- for(int beta=0;beta<2;beta++)
- for(int gamma=0;gamma<2;gamma++)
- for(int delta=0;delta<2;delta++){
- matTvectprod(dta,lb->getphi0(gamma),v1);
- matTvectprod(dtb,lb->getphi0(alpha),v2);
- diaprod(v1,v2,mtemp);
- temp = Hhat->get(alpha,beta,gamma,delta)*(-lb->getphi0(1,delta))*(-lb->getphi0(1,beta));
- for(int i=0;i<3;i++)
- for(int j=0;j<3;j++)
- me[i][j]+=(mtemp[i][j]*temp);
- }
-}
-
-
-inline void BulkC0PlateDGForceBendingTerms(const TensorialTraits<double>::HessType &hessj,const std::vector<TensorialTraits<double>::HessType> &Hess,const LinearElasticShellHookeTensor *H,const LocalBasis *lb,const fullMatrix<double> &disp, double me[3]){
- const int n = Hess.size();
- for(int i=0;i<3;i++) me[i]=0.;
- double mtemp[3][3];
- double B1[3],B2[3];
- for(int a=0;a<2;a++)
- for(int b=0;b<2;b++)
- for(int c=0;c<2;c++)
- for(int d=0;d<2;d++){
- B1[0] = - hessj(c,d)*lb->gett0(0); B1[1] = - hessj(c,d)*lb->gett0(1); B1[2] = - hessj(c,d)*lb->gett0(2);
- for(int j=0;j<n;j++){
- B2[0] = - Hess[j](a,b)*lb->gett0(0); B2[1] = - Hess[j](a,b)*lb->gett0(1); B2[2] = - Hess[j](a,b)*lb->gett0(2);
- diaprod(B1,B2,mtemp);
- for(int i=0;i<3;i++)
- me[i]+=H->get(a,b,c,d)*(mtemp[i][0]*disp(j,0)+mtemp[i][1]*disp(j+n,0)+mtemp[i][2]*disp(j+2*n,0));
- }
- }
-}
-
-inline void BulkC0PlateDGForceMembraneTerms(const double Bj[3][2][2],const double Bn[256][3][2][2] , const LinearElasticShellHookeTensor *H, const int n, const fullMatrix<double> &disp, double me[3]){
- for(int i=0;i<3;i++) me[i]=0.;
- double mtemp[3][3];
- SVector3 v1(0.,0.,0.), v2(0.,0.,0.);
- for(int a=0;a<2;a++)
- for(int b=0;b<2;b++)
- for(int c=0;c<2;c++)
- for(int d=0;d<2;d++){
- v1(0) = Bj[0][c][d]; v1(1) = Bj[1][c][d]; v1[2]=Bj[2][c][d];
- for(int j=0;j<n;j++){
- v2(0) = Bn[j][0][a][b]; v2(1) = Bn[j][1][a][b]; v2(2) = Bn[j][2][a][b];
- diaprod(v1,v2,mtemp);
- for(int jj=0;jj<3;jj++)
- me[jj] += H->get(a,b,c,d)*(mtemp[jj][0]*disp(j,0)+mtemp[jj][1]*disp(j+n,0)+mtemp[jj][2]*disp(j+2*n,0));
- }
- }
-
-}
-
-inline void consC0PlateForceTerms(const LinearElasticShellHookeTensor *Hhat, const int n_m, const int n_p, const double Bhat[3][2][2], const double Dt_m[256][3][3], const double Dt_p[256][3][3],const LocalBasis *lb, const fullMatrix<double> &disp, double me[3]){
- double mtemp[3][3];
- double v1[3],v2[3];
- double temp=0.;
- for(int i=0;i<3;i++) me[i]=0.;
- for(int a=0;a<2;a++)
- for(int b=0;b<2;b++)
- for(int c=0;c<2;c++)
- for(int d=0;d<2;d++){
- v1[0]=Bhat[0][c][d] ; v1[1]=Bhat[1][c][d]; v1[2]=Bhat[2][c][d];
- temp = 0.5*Hhat->get(a,b,c,d)*(-lb->getphi0(1,b));
- for(int j=0;j<n_m;j++){
- matTvectprod(Dt_m[j],lb->getphi0(a),v2); // put in a loop but in this time a vector is need to store the result of matTvectprod
- diaprod(v1,v2,mtemp);
- for(int k=0;k<3;k++)
- me[k] += -temp*(mtemp[k][0]*disp(j,0)+mtemp[k][1]*disp(j+n_m,0)+mtemp[k][2]*disp(j+2*n_m,0));
- }
- for(int j=0;j<n_p;j++){
- matTvectprod(Dt_p[j],lb->getphi0(a),v2); //idem
- diaprod(v1,v2,mtemp);
- for(int k=0;k<3;k++)
- me[k] += temp*(mtemp[k][0]*disp(j+3*n_m,0)+mtemp[k][1]*disp(j+n_p+3*n_m,0)+mtemp[k][2]*disp(j+2*n_p+3*n_m,0));
- }
- }
-}
-
-inline void compC0PlateForceTerms(const int n_m, const int n_p, const LinearElasticShellHookeTensor *Hhat_m, const LinearElasticShellHookeTensor *Hhat_p, const double Bhat_m[256][3][2][2],const double Bhat_p[256][3][2][2],const double Dt[3][3],const LocalBasis *lb, const fullMatrix<double> &disp, double me[3]){
- double v1[3],v2[3];
- double mtemp[3][3];
- double temp_m=0.;
- double temp_p=0.;
- for(int i=0;i<3;i++) me[i]=0.;
- for(int a=0;a<2;a++){
- matTvectprod(Dt,lb->getphi0(a),v1);
- for(int b=0;b<2;b++)
- for(int c=0;c<2;c++)
- for(int d=0;d<2;d++){
- temp_m = 0.5*Hhat_m->get(a,b,c,d)*(-lb->getphi0(1,b));
- temp_p = 0.5*Hhat_p->get(a,b,c,d)*(-lb->getphi0(1,b));
- for(int j=0;j<n_m;j++){
- v2[0] = Bhat_m[j][0][c][d]; v2[1] = Bhat_m[j][1][c][d]; v2[2] = Bhat_m[j][2][c][d];
- diaprod(v1,v2,mtemp);
- for(int jj=0;jj<3;jj++)
- me[jj]+= temp_m*(mtemp[jj][0]*disp(j,0)+mtemp[jj][1]*disp(j+n_m,0)+mtemp[jj][2]*disp(j+2*n_m,0));
- }
- for(int j=0;j<n_p;j++){
- v2[0] = Bhat_p[j][0][c][d]; v2[1] = Bhat_p[j][1][c][d]; v2[2]=Bhat_p[j][2][c][d];
- diaprod(v1,v2,mtemp);
- for(int jj=0;jj<3;jj++)
- me[jj]+= temp_p*(mtemp[jj][0]*disp(j+3*n_m,0)+mtemp[jj][1]*disp(j+n_p+3*n_m,0)+mtemp[jj][2]*disp(j+2*n_p+3*n_m,0));
- }
- }
- }
-}
-
-inline void stabilityC0PlateForceTerms(const int n_p,const int n_m, const LinearElasticShellHookeTensor *Hhat,const double Dt[3][3],const double Dt_m[256][3][3],const double Dt_p[256][3][3], const LocalBasis *lb, const fullMatrix<double> &disp, double me[3]){
- for(int i=0;i<3;i++) me[i]=0.;
- double v1[3],v2[3];
- double mtemp[3][3];
- double temp=0.;
- for(int a=0;a<2;a++)
- for(int b=0;b<2;b++)
- for(int c=0;c<2;c++){
- matTvectprod(Dt,lb->getphi0(c),v1);
- for(int d=0;d<2;d++){
- temp = Hhat->get(a,b,c,d)*(-lb->getphi0(1,d))*(-lb->getphi0(1,b));
- for(int j=0;j<n_m;j++){
- matTvectprod(Dt_m[j],lb->getphi0(a),v2); // can be put on a loop but it is necessary to create a vector to keep the result of matTvectprod
- diaprod(v1,v2,mtemp);
- for(int jj=0;jj<3;jj++)
- me[jj] += -temp*(mtemp[jj][0]*disp(j,0)+mtemp[jj][1]*disp(j+n_m,0)+mtemp[jj][2]*disp(j+2*n_m,0));
- }
- for(int j=0;j<n_p;j++){
- matTvectprod(Dt_p[j],lb->getphi0(a),v2); // idem
- diaprod(v1,v2,mtemp);
- for(int jj=0;jj<3;jj++)
- me[jj] += temp*(mtemp[jj][0]*disp(j+3*n_m,0)+mtemp[jj][1]*disp(j+n_p+3*n_m,0)+mtemp[jj][2]*disp(j+2*n_p+3*n_m,0));
- }
- }
- }
-}
-
-inline void consC0PlateForceTermsBound(const LinearElasticShellHookeTensor *Hhat, const int n, const double Bhat[3][2][2], const double Dt_m[256][3][3], const LocalBasis *lb, const fullMatrix<double> &disp, double me[3]){
- double mtemp[3][3];
- double v1[3],v2[3];
- double temp=0.;
- for(int i=0;i<3;i++) me[i]=0.;
- for(int a=0;a<2;a++)
- for(int b=0;b<2;b++)
- for(int c=0;c<2;c++)
- for(int d=0;d<2;d++){
- v1[0]=Bhat[0][c][d] ; v1[1]=Bhat[1][c][d]; v1[2]=Bhat[2][c][d];
- temp = 0.5*Hhat->get(a,b,c,d)*(-lb->getphi0(1,b));
- for(int j=0;j<n;j++){
- matTvectprod(Dt_m[j],lb->getphi0(a),v2); // put in a loop but in this time a vector is need to store the result of matTvectprod
- diaprod(v1,v2,mtemp);
- for(int k=0;k<3;k++)
- me[k] += -temp*(mtemp[k][0]*disp(j,0)+mtemp[k][1]*disp(j+n,0)+mtemp[k][2]*disp(j+2*n,0));
- }
- }
-}
-
-inline void compC0PlateForceTermsBound(const int n_m,const LinearElasticShellHookeTensor *Hhat_m, const double Bhat_m[256][3][2][2], const double Dt[3][3],const LocalBasis *lb, const fullMatrix<double> &disp, double me[3]){
- double v1[3],v2[3];
- double mtemp[3][3];
- double temp_m=0.;
- for(int i=0;i<3;i++) me[i]=0.;
- for(int a=0;a<2;a++){
- matTvectprod(Dt,lb->getphi0(a),v1);
- for(int b=0;b<2;b++)
- for(int c=0;c<2;c++)
- for(int d=0;d<2;d++){
- temp_m = 0.5*Hhat_m->get(a,b,c,d)*(-lb->getphi0(1,b));
- for(int j=0;j<n_m;j++){
- v2[0] = Bhat_m[j][0][c][d]; v2[1] = Bhat_m[j][1][c][d]; v2[2] = Bhat_m[j][2][c][d];
- diaprod(v1,v2,mtemp);
- for(int jj=0;jj<3;jj++)
- me[jj]+= temp_m*(mtemp[jj][0]*disp(j,0)+mtemp[jj][1]*disp(j+n_m,0)+mtemp[jj][2]*disp(j+2*n_m,0));
- }
- }
- }
-}
-
-inline void stabilityC0PlateForceTermsBound(const int n_m,const LinearElasticShellHookeTensor *Hhat,const double Dt[3][3],const double Dt_m[256][3][3], const LocalBasis *lb, const fullMatrix<double> &disp, double me[3]){
- for(int i=0;i<3;i++) me[i]=0.;
- double v1[3],v2[3];
- double mtemp[3][3];
- double temp=0.;
- for(int a=0;a<2;a++)
- for(int b=0;b<2;b++)
- for(int c=0;c<2;c++){
- matTvectprod(Dt,lb->getphi0(c),v1);
- for(int d=0;d<2;d++){
- temp = Hhat->get(a,b,c,d)*(-lb->getphi0(1,d))*(-lb->getphi0(1,b));
- for(int j=0;j<n_m;j++){
- matTvectprod(Dt_m[j],lb->getphi0(a),v2); // can be put on a loop but it is necessary to create a vector to keep the result of matTvectprod
- diaprod(v1,v2,mtemp);
- for(int jj=0;jj<3;jj++)
- me[jj] += -temp*(mtemp[jj][0]*disp(j,0)+mtemp[jj][1]*disp(j+n_m,0)+mtemp[jj][2]*disp(j+2*n_m,0));
- }
- }
- }
-}
-
-// Compute value needed in get
-void compute_Bn(const LocalBasis *lb, const std::vector<TensorialTraits<SVector3>::GradType> &Grads, const int n, double B[][3][2][2]){
- for(int i=0;i<n;i++){
- for(int a=0;a<2;a++)
- for(int b=0;b<2;b++){
- B[i][0][a][b] = 0.5*(lb->getphi0(a,0)*Grads[i](0,b)+lb->getphi0(b,0)*Grads[i](0,a) );
- B[i][1][a][b] = 0.5*(lb->getphi0(a,1)*Grads[i](0,b)+lb->getphi0(b,1)*Grads[i](0,a) );
- B[i][2][a][b] = 0.5*(lb->getphi0(a,2)*Grads[i](0,b)+lb->getphi0(b,2)*Grads[i](0,a) );
- }
- }
-
-}
-
-void Compute_Bhat(const LocalBasis *lb,const std::vector<TensorialTraits<double>::HessType> &Hess, const int &n, double B[256][3][2][2]){
- for(int i=0;i<n;i++){
- for(int ii=0;ii<3;ii++){
- B[i][ii][0][0] =0.; B[i][ii][0][1]=0.; B[i][ii][1][0]=0.; B[i][ii][1][1]=0.;
- for(int j=0;j<2;j++)
- for(int k=0;k<2;k++){
- B[i][ii][0][0] += -lb->gett0(ii)*lb->gett(0,j)*Hess[i](j,k)*lb->gett(0,k);
- B[i][ii][0][1] += -lb->gett0(ii)*lb->gett(0,j)*Hess[i](j,k)*lb->gett(1,k);
- B[i][ii][1][0] += -lb->gett0(ii)*lb->gett(1,j)*Hess[i](j,k)*lb->gett(0,k);
- B[i][ii][1][1] += -lb->gett0(ii)*lb->gett(1,j)*Hess[i](j,k)*lb->gett(1,k);
- }
- }
- }
-}
-
-// Compute value needed in getinter
-void compute_Deltat_tilde(const LocalBasis *lb, const std::vector<TensorialTraits<SVector3>::GradType> &Grads, const int &n, double Deltat[256][3][3]){
- SVector3 vect(0.,0.,0.),vect2(0.,0.,0.);
- vect = crossprod(lb->gett0(),lb->getphi0(0));
- vect2= crossprod(lb->gett0(),lb->getphi0(1));
- double invJ0 = 1./lb->getJacobian();
- for(int i=0;i<n;i++){
- Deltat[i][0][0] = invJ0 *(-lb->gett0(0)*vect(2)*Grads[i](0,1) + lb->gett0(0)*vect2(2))*Grads[i](0,0);
- Deltat[i][1][0] = invJ0 *((lb->getphi0(0,2)-lb->gett0(1)*vect(2))*Grads[i](0,1) - (lb->getphi0(1,2)-lb->gett0(1)*vect2(2))*Grads[i](0,0));
- Deltat[i][2][0] = invJ0 *((-lb->getphi0(0,1)-lb->gett0(2)*vect(2))*Grads[i](0,1) - (-lb->getphi0(1,1)-lb->gett0(2)*vect2(2))*Grads[i](0,0));
-
- Deltat[i][0][1] = invJ0 *((-lb->getphi0(0,2)-lb->gett0(0)*vect(2))*Grads[i](0,1) - (-lb->getphi0(1,2)-lb->gett0(0)*vect2(2))*Grads[i](0,0));
- Deltat[i][1][1] = invJ0 *(-lb->gett0(1)*vect(2)*Grads[i](0,1) + lb->gett0(1)*vect2(2)*Grads[i](0,0));
- Deltat[i][2][1] = invJ0 *((lb->getphi0(0,0)-lb->gett0(2)*vect(2))*Grads[i](0,1) - (lb->getphi0(1,0)-lb->gett0(2)*vect2(2))*Grads[i](0,0));
-
- Deltat[i][0][2] = invJ0 *((lb->getphi0(0,1)-lb->gett0(0)*vect(2))*Grads[i](0,1) - (lb->getphi0(1,1)-lb->gett0(0)*vect2(2))*Grads[i](0,0));
- Deltat[i][1][2] = invJ0 *((-lb->getphi0(0,0)-lb->gett0(1)*vect(2))*Grads[i](0,1) - (-lb->getphi0(1,0)-lb->gett0(1)*vect2(2))*Grads[i](0,0));
- Deltat[i][2][2] = invJ0 *( -lb->gett0(2)*vect(2)*Grads[i](0,1) + lb->gett0(2)*vect2(2)*Grads[i](0,0));
- }
-}
-void compute_Deltat_tildeBound(const LocalBasis *lb, const std::vector<TensorialTraits<SVector3>::GradType> &Grads, const int &n, double Deltat[256][3][3], const LocalBasis *lbs){
- SVector3 vect(0.,0.,0.),vect2(0.,0.,0.),temp2(0.,0.,0.),temp3(0.,0.,0.),phi1norm(0.,0.,0.);
- vect = crossprod(lb->gett0(),lb->getphi0(0));
- vect2= crossprod(lb->gett0(),lb->getphi0(1));
- phi1norm=lbs->getphi0(0);
- phi1norm.normalize();
- double temp[3][3];
- double invJ0 = 1./lb->getJacobian();
- for(int i=0;i<n;i++){
- temp[0][0] = invJ0 *(-lb->gett0(0)*vect(2)*Grads[i](0,1) + lb->gett0(0)*vect2(2))*Grads[i](0,0);
- temp[1][0] = invJ0 *((lb->getphi0(0,2)-lb->gett0(1)*vect(2))*Grads[i](0,1) - (lb->getphi0(1,2)-lb->gett0(1)*vect2(2))*Grads[i](0,0));
- temp[2][0] = invJ0 *((-lb->getphi0(0,1)-lb->gett0(2)*vect(2))*Grads[i](0,1) - (-lb->getphi0(1,1)-lb->gett0(2)*vect2(2))*Grads[i](0,0));
-
- temp[0][1] = invJ0 *((-lb->getphi0(0,2)-lb->gett0(0)*vect(2))*Grads[i](0,1) - (-lb->getphi0(1,2)-lb->gett0(0)*vect2(2))*Grads[i](0,0));
- temp[1][1] = invJ0 *(-lb->gett0(1)*vect(2)*Grads[i](0,1) + lb->gett0(1)*vect2(2)*Grads[i](0,0));
- temp[2][1] = invJ0 *((lb->getphi0(0,0)-lb->gett0(2)*vect(2))*Grads[i](0,1) - (lb->getphi0(1,0)-lb->gett0(2)*vect2(2))*Grads[i](0,0));
-
- temp[0][2] = invJ0 *((lb->getphi0(0,1)-lb->gett0(0)*vect(2))*Grads[i](0,1) - (lb->getphi0(1,1)-lb->gett0(0)*vect2(2))*Grads[i](0,0));
- temp[1][2] = invJ0 *((-lb->getphi0(0,0)-lb->gett0(1)*vect(2))*Grads[i](0,1) - (-lb->getphi0(1,0)-lb->gett0(1)*vect2(2))*Grads[i](0,0));
- temp[2][2] = invJ0 *( -lb->gett0(2)*vect(2)*Grads[i](0,1) + lb->gett0(2)*vect2(2)*Grads[i](0,0));
-
- // Projection of normal
- for(int j=0;j<3;j++){
- temp3(0)=temp[0][j]; temp3(1)=temp[1][j]; temp3(2)=temp[2][j];
- temp2=dot(temp3,phi1norm)*phi1norm;
- temp3-=temp2;
- Deltat[i][0][j]=temp3(0);Deltat[i][1][j]=temp3(1); Deltat[i][2][j]=temp3(2);
- }
- }
-}
-#endif //DGC0PLATEELEMENTARYTERMS_H_
diff --git a/contrib/cm3/DgC0PlateSolver.cpp b/contrib/cm3/DgC0PlateSolver.cpp
deleted file mode 100644
index b4bbd36a29..0000000000
--- a/contrib/cm3/DgC0PlateSolver.cpp
+++ /dev/null
@@ -1,386 +0,0 @@
-// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to <gmsh@geuz.org>.
-
-#include <string.h>
-#include "GmshConfig.h"
-#include "DgC0PlateSolver.h"
-#include "linearSystemCSR.h"
-#include "linearSystemPETSc.h"
-#include "linearSystemGMM.h"
-#include "Numeric.h"
-#include "GModelWithInterface.h"
-#include "functionSpace.h"
-#include "C0DgPlateTerms.h"
-#include "solverAlgorithms.h"
-#include "quadratureRules.h"
-#include "solverField.h"
-#include "AssembleInterface.h"
-#include "MPoint.h"
-
-#if defined(HAVE_POST)
-#include "PView.h"
-#include "PViewData.h"
-#endif
-
-void DgC0PlateSolver::setMesh(const std::string &meshFileName)
-{
- pModel = new GModelWithInterface();
- pModel->readMSH(meshFileName.c_str());
- _dim = pModel->getNumRegions() ? 3 : 2;
- if (LagSpace) delete LagSpace;
- // Faudra quelque chose dans le fichier de données qui permettra de choisir
- //if (_dim==3) LagSpace=new VectorLagrangeFunctionSpace(_tag);
- //if (_dim==2) LagSpace=new VectorLagrangeFunctionSpace(_tag,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y);
- // Plaque un ddl par noeud
- //LagSpace=new VectorLagrangeFunctionSpace(_tag,VectorLagrangeFunctionSpace::VECTOR_Z);
- // Plate 3 dof per node
- LagSpace=new VectorLagrangeFunctionSpace(_tag);
-}
-
-void DgC0PlateSolver::readInputFile(const std::string &fn)
-{
- FILE *f = fopen(fn.c_str(), "r");
- char what[256];
- while(!feof(f)){
- if(fscanf(f, "%s", what) != 1) return;
- if (!strcmp(what, "ElasticDomain")){
- DGelasticField field;
- int physical;
- // Add two parameters beta1 (used for stabilization) and h the height of the plate
- if(fscanf(f, "%d %lf %lf %lf %lf", &physical, &field._E, &field._nu, &field._beta1, &field._h) != 5) return;
- field._tag = _tag;
- field.g = new groupOfElements (_dim, physical);
- // Add the Interface Elements
- field.gi = pModel->getInterface(physical); // TODO integrate with field.g
- elasticFields.push_back(field);
- }
- else if (!strcmp(what, "Void")){
- // elasticField field;
- // create enrichment ...
- // create the group ...
- // assign a tag
- // elasticFields.push_back(field);
- }
- else if (!strcmp(what, "NodeDisplacement")){
- double val;
- int node, comp;
- if(fscanf(f, "%d %d %lf", &node, &comp, &val) != 3) return;
- dirichletBC diri;
- diri.g = new groupOfElements (0, node);
- diri._f= simpleFunction<double>(val);
- diri._comp=comp;
- diri._tag=node;
- diri.onWhat=BoundaryCondition::ON_VERTEX;
- allDirichlet.push_back(diri);
- }
- else if (!strcmp(what, "EdgeDisplacement")){
- double val;
- int edge, comp;
- if(fscanf(f, "%d %d %lf", &edge, &comp, &val) != 3) return;
- dirichletBC diri;
- diri.g = new groupOfElements (1, edge);
- diri._f= simpleFunction<double>(val);
- diri._comp=comp;
- diri._tag=edge;
- diri.onWhat=BoundaryCondition::ON_EDGE;
- allDirichlet.push_back(diri);
- }
- else if (!strcmp(what, "FaceDisplacement")){
- double val;
- int face, comp;
- if(fscanf(f, "%d %d %lf", &face, &comp, &val) != 3) return;
- dirichletBC diri;
- diri.g = new groupOfElements (2, face);
- diri._f= simpleFunction<double>(val);
- diri._comp=comp;
- diri._tag=face;
- diri.onWhat=BoundaryCondition::ON_FACE;
- allDirichlet.push_back(diri);
- }
- else if (!strcmp(what, "NodeForce")){
- double val1, val2, val3;
- int node;
- if(fscanf(f, "%d %lf %lf %lf", &node, &val1, &val2, &val3) != 4) return;
- neumannBC neu;
- neu.g = new groupOfElements (0, node);
- neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
- neu._tag=node;
- neu.onWhat=BoundaryCondition::ON_VERTEX;
- allNeumann.push_back(neu);
- }
- else if (!strcmp(what, "EdgeForce")){
- double val1, val2, val3;
- int edge;
- if(fscanf(f, "%d %lf %lf %lf", &edge, &val1, &val2, &val3) != 4) return;
- neumannBC neu;
- neu.g = new groupOfElements (1, edge);
- neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
- neu._tag=edge;
- neu.onWhat=BoundaryCondition::ON_EDGE;
- allNeumann.push_back(neu);
- }
- else if (!strcmp(what, "FaceForce")){
- double val1, val2, val3;
- int face;
- if(fscanf(f, "%d %lf %lf %lf", &face, &val1, &val2, &val3) != 4) return;
- neumannBC neu;
- neu.g = new groupOfElements (2, face);
- neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
- neu._tag=face;
- neu.onWhat=BoundaryCondition::ON_FACE;
- allNeumann.push_back(neu);
- }
- else if (!strcmp(what, "VolumeForce")){
- double val1, val2, val3;
- int volume;
- if(fscanf(f, "%d %lf %lf %lf", &volume, &val1, &val2, &val3) != 4) return;
- neumannBC neu;
- neu.g = new groupOfElements (3, volume);
- neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
- neu._tag=volume;
- neu.onWhat=BoundaryCondition::ON_VOLUME;
- allNeumann.push_back(neu);
- }
- else if (!strcmp(what, "MeshFile")){
- char name[245];
- if(fscanf(f, "%s", name) != 1) return;
- setMesh(name);
- }
- else if (!strcmp(what, "Theta")){
- int num, num_phys=0;
- // First component number of physical groups where theta is imposed (to 0 for now)
- fscanf(f,"%d",&num);
- for(int i=0;i<num;i++){
- fscanf(f,"%d",&num_phys);
- // Find the interface element linked to the physical group num_phys (one interfaceelement vector for all physical)
- pModel->generateInterfaceElementsOnBoundary(num_phys,elasticFields); //TODO don't pass elastic field but I don't know how to access to element in GModel ??
- }
- // store the boundary interface element to the elastifFields
- for(int i=0;i<elasticFields.size();i++) elasticFields[i].gib=pModel->getBoundInterface(i);
- }
- else {
- Msg::Error("Invalid input : %s", what);
-// return;
- }
- }
- fclose(f);
-}
-
-void DgC0PlateSolver::solve()
-{
-#if defined(HAVE_TAUCS)
- linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
-#elif defined(HAVE_PETSC)
- linearSystemPETSc<double> *lsys = new linearSystemPETSc<double>;
-#else
- linearSystemGmm<double> *lsys = new linearSystemGmm<double>;
- lsys->setNoisy(2);
-#endif
- if (pAssembler) delete pAssembler;
- pAssembler = new dofManager<double>(lsys);
-
- // we first do all fixations. the behavior of the dofManager is to
- // give priority to fixations : when a dof is fixed, it cannot be
- // numbered afterwards
-
- std::cout << "Dirichlet BC"<< std::endl;
- for (unsigned int i = 0; i < allDirichlet.size(); i++)
- {
- FilterDofComponent filter(allDirichlet[i]._comp);
- FixNodalDofs(*LagSpace,allDirichlet[i].g->begin(),allDirichlet[i].g->end(),*pAssembler,allDirichlet[i]._f,filter);
- }
-
- // we number the dofs : when a dof is numbered, it cannot be numbered
- // again with another number.
- for (unsigned int i = 0; i < elasticFields.size(); ++i)
- {
- NumberDofs(*LagSpace, elasticFields[i].g->begin(), elasticFields[i].g->end(),*pAssembler);
- }
-
- // Now we start the assembly process
- // First build the force vector
-
- GaussQuadrature Integ_Boundary(GaussQuadrature::Val);
- std::cout << "Neumann BC"<< std::endl;
- for (unsigned int i = 0; i < allNeumann.size(); i++)
- {
- LoadTerm<SVector3> Lterm(*LagSpace,allNeumann[i]._f);
- Assemble(Lterm,*LagSpace,allNeumann[i].g->begin(),allNeumann[i].g->end(),Integ_Boundary,*pAssembler);
- }
-
-// bulk material law
-
- GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
- for (unsigned int i = 0; i < elasticFields.size(); i++)
- {
- // Initialization of elementary terms in function of the field and space
- IsotropicElasticTermC0Plate Eterm(*LagSpace,elasticFields[i]._E,elasticFields[i]._nu,elasticFields[i]._h);
-
- // Assembling loop on Elementary terms
- Assemble(Eterm,*LagSpace,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,*pAssembler);
-
- // Initialization of elementary interface terms in function of the field and space
- IsotropicElasticInterfaceTermC0Plate IEterm(*LagSpace,elasticFields[i]._E,elasticFields[i]._nu,elasticFields[i]._beta1,elasticFields[i]._h);
-
- // Assembling loop on elementary interface terms
- AssembleInterface(IEterm,*LagSpace,elasticFields[i].g,elasticFields[i].gi,Integ_Boundary,*pAssembler); // Use the same GaussQuadrature rule than on the boundary
-
- // Assembling loop on elementary boundary interface terms
- AssembleInterface(IEterm,*LagSpace,elasticFields[i].g,elasticFields[i].gib,Integ_Boundary,*pAssembler); // Use the same GaussQuadrature rule than on the boundary
-
- }
-printf("-- done assembling!\n");
-lsys->systemSolve();
-printf("-- done solving!\n");
-
- double energ=0;
- for (unsigned int i = 0; i < elasticFields.size(); i++)
- {
- SolverField<SVector3> Field(pAssembler, LagSpace);
- IsotropicElasticTermC0Plate Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu,elasticFields[i]._h);
- BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
- Assemble(Elastic_Energy_Term,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,energ);
- }
- printf("elastic energy=%f\n",energ);
-}
-
-
-
-#if defined(HAVE_POST)
-
-static double vonMises(dofManager<double> *a, MElement *e,
- double u, double v, double w,
- double E, double nu, int _tag)
-{
- double valx[256];
- double valy[256];
- double valz[256];
- for (int k = 0; k < e->getNumVertices(); k++){
- valx[k] = a->getDofValue(e->getVertex(k), 0, _tag);
- valy[k] = a->getDofValue(e->getVertex(k), 1, _tag);
- valz[k] = a->getDofValue(e->getVertex(k), 2, _tag);
- }
- double gradux[3];
- double graduy[3];
- double graduz[3];
- e->interpolateGrad(valx, u, v, w, gradux);
- e->interpolateGrad(valy, u, v, w, graduy);
- e->interpolateGrad(valz, u, v, w, graduz);
-
- double eps[6] = {gradux[0], graduy[1], graduz[2],
- 0.5 * (gradux[1] + graduy[0]),
- 0.5 * (gradux[2] + graduz[0]),
- 0.5 * (graduy[2] + graduz[1])};
- double A = E / (1. + nu);
- double B = A * (nu / (1. - 2 * nu));
- double trace = eps[0] + eps[1] + eps[2] ;
- double sxx = A * eps[0] + B * trace;
- double syy = A * eps[1] + B * trace;
- double szz = A * eps[2] + B * trace;
- double sxy = A * eps[3];
- double sxz = A * eps[4];
- double syz = A * eps[5];
-
- double s[9] = {sxx, sxy, sxz,
- sxy, syy, syz,
- sxz, syz, szz};
-
- return ComputeVonMises(s);
-}
-
-PView* DgC0PlateSolver::buildDisplacementView (const std::string &postFileName)
-{
- std::set<MVertex*> v;
- for (unsigned int i = 0; i < elasticFields.size(); ++i)
- {
- for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it)
- {
- MElement *e=*it;
- for (int j = 0; j < e->getNumVertices(); ++j) v.insert(e->getVertex(j));
- }
- }
- std::map<int, std::vector<double> > data;
- SolverField<SVector3> Field(pAssembler, LagSpace);
- for ( std::set<MVertex*>::iterator it = v.begin(); it != v.end(); ++it)
- {
- SVector3 val;
- MPoint p(*it);
- Field.f(&p,0,0,0,val);
- std::vector<double> vec(3);vec[0]=val(0);vec[1]=val(1);vec[2]=val(2);
- data[(*it)->getNum()]=vec;
- }
- PView *pv = new PView (postFileName, "NodeData", pModel, data, 0.0);
- return pv;
-}
-
-PView *DgC0PlateSolver::buildElasticEnergyView(const std::string &postFileName)
-{
- std::map<int, std::vector<double> > data;
- GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
- for (unsigned int i = 0; i < elasticFields.size(); ++i)
- {
- SolverField<SVector3> Field(pAssembler, LagSpace);
- IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
- BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
- ScalarTermConstant One(1.0);
- for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it)
- {
- MElement *e=*it;
- double energ;
- double vol;
- IntPt *GP;
- int npts=Integ_Bulk.getIntPoints(e,&GP);
- Elastic_Energy_Term.get(e,npts,GP,energ);
- One.get(e,npts,GP,vol);
- std::vector<double> vec;
- vec.push_back(energ/vol);
- data[(*it)->getNum()]=vec;
- }
- }
- PView *pv = new PView (postFileName, "ElementData", pModel, data, 0.0);
- return pv;
-}
-
-PView *DgC0PlateSolver::buildVonMisesView(const std::string &postFileName)
-{
- std::map<int, std::vector<double> > data;
- GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
- for (unsigned int i = 0; i < elasticFields.size(); ++i)
- {
- SolverField<SVector3> Field(pAssembler, LagSpace);
- IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
- BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
- for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it)
- {
- MElement *e=*it;
- double energ;
- double vol;
- IntPt *GP;
- int npts=Integ_Bulk.getIntPoints(e,&GP);
- Elastic_Energy_Term.get(e,npts,GP,energ);
- std::vector<double> vec;
- vec.push_back(energ);
- data[(*it)->getNum()]=vec;
- }
- }
- PView *pv = new PView (postFileName, "ElementData", pModel, data, 0.0);
- return pv;
-}
-
-
-#else
-PView* DgC0PlateSolver::buildDisplacementView (const std::string &postFileName)
-{
- Msg::Error("Post-pro module not available");
- return 0;
-}
-
-PView* DgC0PlateSolver::buildElasticEnergyView(const std::string &postFileName)
-{
- Msg::Error("Post-pro module not available");
- return 0;
-}
-
-#endif
diff --git a/contrib/cm3/DgC0PlateSolver.h b/contrib/cm3/DgC0PlateSolver.h
deleted file mode 100644
index 0d3689426b..0000000000
--- a/contrib/cm3/DgC0PlateSolver.h
+++ /dev/null
@@ -1,96 +0,0 @@
-// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to <gmsh@geuz.org>.
-
-#ifndef _DGC0PLATE_SOLVER_H_
-#define _DGC0PLATE_SOLVER_H_
-
-#include <map>
-#include <string>
-#include "SVector3.h"
-#include "dofManager.h"
-#include "simpleFunction.h"
-#include "functionSpace.h"
-#include "MInterfaceElement.h"
-
-class GModelWithInterface;
-class PView;
-class groupOfElements;
-
-struct elasticField {
- int _tag; // tag for the dofManager
- groupOfElements *g; // support for this field
- double _E, _nu; // specific elastic datas (should be somewhere else)
- elasticField () : g(0),_tag(0){}
-};
-
-// Elasticfield for DG element (has a group of elements for interfaceelement)
-struct DGelasticField{
- int _tag; // tag for the dofManager
- groupOfElements *g; // support for this field
- std::vector<MInterfaceElement*> gi; // support for the interfaceElement TODO cast to a groupOfElements
- std::vector<MInterfaceElement*> gib; // support for the interfaceElement TODO cast to a groupOfElements
- double _E, _nu, _beta1, _h; // specific elastic datas (should be somewhere else)
- DGelasticField () : g(0), _tag(0){}
-};
-
-struct BoundaryCondition
-{
- enum location{UNDEF,ON_VERTEX,ON_EDGE,ON_FACE,ON_VOLUME};
- location onWhat; // on vertices or elements
- int _tag; // tag for the dofManager
- groupOfElements *g; // support for this BC
- BoundaryCondition() : g(0),_tag(0),onWhat(UNDEF) {};
-};
-
-struct dirichletBC : public BoundaryCondition
-{
- int _comp; // component
- simpleFunction<double> _f;
- dirichletBC () :BoundaryCondition(),_comp(0),_f(0){}
-};
-
-struct neumannBC : public BoundaryCondition
-{
- simpleFunction<SVector3> _f;
- neumannBC () : BoundaryCondition(),_f(SVector3(0,0,0)){}
-};
-
-// an elastic solver ...
-class DgC0PlateSolver
-{
- protected:
- GModelWithInterface *pModel;
- int _dim, _tag;
- dofManager<double> *pAssembler;
- FunctionSpace<SVector3> *LagSpace;
-
- // young modulus and poisson coefficient per physical
- std::vector<DGelasticField> elasticFields;
- // neumann BC
- std::vector<neumannBC> allNeumann;
- // dirichlet BC
- std::vector<dirichletBC> allDirichlet;
-
- public:
- DgC0PlateSolver(int tag) : _tag(tag),LagSpace(0),pAssembler(0) {}
- virtual ~DgC0PlateSolver()
- {
- if (LagSpace) delete LagSpace;
- if (pAssembler) delete pAssembler;
- }
- void readInputFile(const std::string &meshFileName);
- void setMesh(const std::string &meshFileName);
- virtual void solve();
- virtual PView *buildDisplacementView(const std::string &postFileName);
- virtual PView *buildElasticEnergyView(const std::string &postFileName);
- virtual PView *buildVonMisesView(const std::string &postFileName);
- // std::pair<PView *, PView*> buildErrorEstimateView
- // (const std::string &errorFileName, double, int);
- // std::pair<PView *, PView*> buildErrorEstimateView
- // (const std::string &errorFileName, const elasticityData &ref, double, int);
-};
-
-
-#endif
diff --git a/contrib/cm3/GModelWithInterface.cpp b/contrib/cm3/GModelWithInterface.cpp
deleted file mode 100644
index cad85d7349..0000000000
--- a/contrib/cm3/GModelWithInterface.cpp
+++ /dev/null
@@ -1,664 +0,0 @@
-//
-// C++ Interface: terms
-//
-// Description: Function of GModelWithInterface
-//
-//
-// Author: <Gauthier BECKER>, (C) 2010
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-
-#include "GModelWithInterface.h"
-#include "GEntity.h"
-#include "MInterfaceElement.h"
-#include <limits>
-#include <stdlib.h>
-#include <string.h>
-#include <map>
-#include <string>
-#include "GmshDefines.h"
-#include "MPoint.h"
-#include "MLine.h"
-#include "MTriangle.h"
-#include "MQuadrangle.h"
-#include "MTetrahedron.h"
-#include "MHexahedron.h"
-#include "MPrism.h"
-#include "MPyramid.h"
-#include "MElementCut.h"
-#include "SBoundingBox3d.h"
-#include "StringUtils.h"
-#include "GmshMessage.h"
-#include "discreteVertex.h"
-#include "discreteEdge.h"
-#include "discreteFace.h"
-#include "discreteRegion.h"
-#include "MElement.h"
-#include "GEdgeCompound.h"
-#include "GFaceCompound.h"
-
-
-// Function allowing to create an interface element TODO try to optimized it (use the work done for interelement on boundary)
-// TODO verification of this function with element plus of polynomialorder more important than minus element
-void createInterfaceElementMSH(GModel *m, MElement *e,
- int reg, int part, std::vector<MVertex*> &v,
- std::map<int, std::vector<MElement*> > elements[10],
- std::vector<std::pair<int,MInterfaceElement> > &interelements, const int physicals)
-{
-
- // take the vector of element corresponding to the element's type
- std::vector<MElement*> ve;
- switch(e->getType()){
- case TYPE_PNT :
- ve = elements[0][reg]; break;
- case TYPE_LIN :
- ve = elements[1][reg]; break;
- case TYPE_TRI :
- ve = elements[2][reg]; break;
- case TYPE_QUA :
- ve = elements[3][reg]; break;
- case TYPE_TET :
- ve = elements[4][reg]; break;
- case TYPE_HEX :
- ve = elements[5][reg]; break;
- case TYPE_PRI :
- ve = elements[6][reg]; break;
- case TYPE_PYR :
- ve = elements[7][reg]; break;
- case TYPE_POLYG :
- ve = elements[8][reg]; break;
- case TYPE_POLYH :
- ve = elements[9][reg]; break;
- default : Msg::Error("Wrong type of element");
- }
-
- // (Primary) Nodes of element e
- int npn = e->getNumPrimaryVertices();
- std::vector<MVertex*> nodes;
-// std::vector<MVertex*> *nodes = NULL;
- for (int i = 0; i < npn ; i++)
- nodes.push_back(e->getVertex(i));
- // loop on element to find the interface elements associate to this element (begin at the end of the vector because the elements
- // adjoining the element are a number near the looked element )
- MElement* el = NULL;
- std::vector<MVertex*> nodes_el;
- MVertex* pt;
- double dist = 0.;
- int count = 0; // count the number of interface elements found
- for( int i = ve.size() - 2; i > -1; i--) // -2 because the last element has no interface element with itself
- {
- el=ve[i];
- // List of el's (primary) vertices
- nodes_el.clear();
- for(int j = 0; j < el->getNumPrimaryVertices(); j++)
- nodes_el.push_back(el->getVertex(j));
- int first_point_m = -1; // First point of interface element has number of vertex is >=0 it can be initialized to -1
- int first_point_p = -1;
-
- // loop on el's (primary) vertices // TODO The last node is not considered because if the others is negative the last one will be negative
- for(int j = 0; j<nodes_el.size(); j++){
- pt = nodes_el[j];
-
- for(int k = 0; k < nodes.size(); k++){
- if (pt == nodes[k]){
- if (first_point_m == -1){
- first_point_m = k;
- first_point_p = j;
- break;
- }
- else{
- // Creation of interface element
- // It is a MLine of the degree of an edge of the element --> the vertices of the (minus) element must be retrieve
- std::vector<MVertex*> vv;
- // choose if the node of plus and minus element that are kept
- int diffpo = e->getPolynomialOrder() - el->getPolynomialOrder();
- if(diffpo <0){
- switch(el->getType()){
- case TYPE_TRI :
- switch(first_point_p+j){
- case 1 :
- el->getEdgeVertices(0,vv); break;
- case 2 :
- el->getEdgeVertices(2,vv); break;
- case 3 :
- el->getEdgeVertices(1,vv); break;
- }
- break;
- case TYPE_QUA :
- switch(first_point_p+j){
- case 1 :
- el->getEdgeVertices(0,vv); break;
- case 5 :
- el->getEdgeVertices(2,vv); break;
- case 3 :
- if(first_point_p==0 or k==0) el->getEdgeVertices(3,vv);
- else el->getEdgeVertices(2,vv);
- break;
- }
- break;
- default : Msg::Error("Cannot build interface element for element type : %s",e->getType());
- }
- MInterfaceElement ie = MInterfaceElement(vv, 0, part, el, e); // MinterfaceElement defined with the same direction of minus element
- interelements.push_back(std::pair<int,MInterfaceElement>(physicals,ie));
- }
- else{
- switch(e->getType()){
- case TYPE_TRI :
- switch(first_point_m+k){
- case 1 :
- e->getEdgeVertices(0,vv); break;
- case 2 :
- e->getEdgeVertices(2,vv); break;
- case 3 :
- e->getEdgeVertices(1,vv); break;
- }
- break;
- case TYPE_QUA :
- switch(first_point_m+k){
- case 1 :
- e->getEdgeVertices(0,vv); break;
- case 5 :
- e->getEdgeVertices(2,vv); break;
- case 3 :
- if(first_point_m==0 or k==0) e->getEdgeVertices(3,vv);
- else e->getEdgeVertices(2,vv);
- break;
- }
- break;
- default : Msg::Error("Cannot build interface element for element type : %s",e->getType());
- }
- MInterfaceElement ie = MInterfaceElement(vv, 0, part, e, el); // MinterfaceElement defined with the same direction of minus element
- interelements.push_back(std::pair<int,MInterfaceElement>(physicals,ie));
- }
- vv.clear();
- count++;
- break;
- }
- }
- }
- }
- // Check if the maximum interface element is reached
- if(count == npn) break;
- }
-}
-
-
-// Function to store the interfaceElement in the GModel
-//void GModelWithInterface::_storeInterfaceElements(std::vector<std::map<int,MInterfaceElement> > & ie )
-//{
-// interfaces=ie;
-//}
-
-
-
-// Duplication of GModelIO_Mesh.cpp because it's static functions ??
-static bool getVertices(int num, int *indices, std::map<int, MVertex*> &map,
- std::vector<MVertex*> &vertices)
-{
- for(int i = 0; i < num; i++){
- if(!map.count(indices[i])){
- Msg::Error("Wrong vertex index %d", indices[i]);
- return false;
- }
- else
- vertices.push_back(map[indices[i]]);
- }
- return true;
-}
-
-static bool getVertices(int num, int *indices, std::vector<MVertex*> &vec,
- std::vector<MVertex*> &vertices)
-{
- for(int i = 0; i < num; i++){
- if(indices[i] < 0 || indices[i] > (int)(vec.size() - 1)){
- Msg::Error("Wrong vertex index %d", indices[i]);
- return false;
- }
- else
- vertices.push_back(vec[indices[i]]);
- }
- return true;
-}
-
-static MElement *createElementMSH(GModel *m, int num, int typeMSH, int physical,
- int reg, int part, std::vector<MVertex*> &v,
- std::map<int, std::vector<MElement*> > elements[10],
- std::map<int, std::map<int, std::string> > physicals[4])
-{
- MElementFactory factory;
- MElement *e = factory.create(typeMSH, v, num, part);
-
- if(!e){
- Msg::Error("Unknown type of element %d", typeMSH);
- return NULL;
- }
-
- switch(e->getType()){
- case TYPE_PNT :
- elements[0][reg].push_back(e); break;
- case TYPE_LIN :
- elements[1][reg].push_back(e); break;
- case TYPE_TRI :
- elements[2][reg].push_back(e); break;
- case TYPE_QUA :
- elements[3][reg].push_back(e); break;
- case TYPE_TET :
- elements[4][reg].push_back(e); break;
- case TYPE_HEX :
- elements[5][reg].push_back(e); break;
- case TYPE_PRI :
- elements[6][reg].push_back(e); break;
- case TYPE_PYR :
- elements[7][reg].push_back(e); break;
- case TYPE_POLYG :
- elements[8][reg].push_back(e); break;
- case TYPE_POLYH :
- elements[9][reg].push_back(e); break;
- default : Msg::Error("Wrong type of element"); return NULL;
- }
-
- int dim = e->getDim();
- if(physical && (!physicals[dim].count(reg) || !physicals[dim][reg].count(physical)))
- physicals[dim][reg][physical] = "unnamed";
-
- if(part) m->getMeshPartitions().insert(part);
- return e;
-}
-
-static MElement *getParent(int parentNum, std::map<int, std::vector<MElement*> > &elements)
-{
- std::map<int, std::vector<MElement*> >::iterator it = elements.begin();
- for(; it != elements.end(); ++it) {
- std::vector<MElement*>::iterator itE = it->second.begin();
- for(; itE != it->second.end(); itE++)
- if((*itE)->getNum() == parentNum) {
- MElement *p = (*itE);
- it->second.erase(itE);
- return p;
- }
- }
- return NULL;
-}
-
-static MElement *getParent(int parentNum, int dim,
- std::map<int, std::vector<MElement*> > elements[10])
-{
- MElement *parent = NULL;
- switch(dim){
- case 0 :
- return getParent(parentNum, elements[0]);
- case 1 :
- return getParent(parentNum, elements[1]);
- case 2 :
- if(parent = getParent(parentNum, elements[2])) return parent;
- if(parent = getParent(parentNum, elements[3])) return parent;
- if(parent = getParent(parentNum, elements[8])) return parent;
- return parent;
- case 3 :
- if(parent = getParent(parentNum, elements[4])) return parent;
- if(parent = getParent(parentNum, elements[5])) return parent;
- if(parent = getParent(parentNum, elements[6])) return parent;
- if(parent = getParent(parentNum, elements[7])) return parent;
- if(parent = getParent(parentNum, elements[9])) return parent;
- return parent;
- default : return parent;
- }
-}
-
-
-// Read the mesh from .msh file +/- a copy of function of GModel but a line is added to create interfaceelemnt and a if(THETA) is added to compute interface element on extremities
-int GModelWithInterface::readMSH(const std::string &name)
-{
- FILE *fp = fopen(name.c_str(), "rb");
- if(!fp){
- Msg::Error("Unable to open file '%s'", name.c_str());
- return 0;
- }
-
- char str[256] = "XXX";
- double version = 1.0;
- bool binary = false, swap = false, postpro = false;
- std::map<int, std::vector<MElement*> > elements[10];
- //std::vector<std::map<int,MInterfaceElement> > interelements;
- std::map<int, std::map<int, std::string> > physicals[4];
- std::map<int, MVertex*> vertexMap;
- std::vector<MVertex*> vertexVector;
-
-
- while(1) {
-
- while(str[0] != '$'){
- if(!fgets(str, sizeof(str), fp) || feof(fp))
- break;
- }
-
- if(feof(fp))
- break;
-
- if(!strncmp(&str[1], "MeshFormat", 10)) {
-
- if(!fgets(str, sizeof(str), fp)) return 0;
- int format, size;
- if(sscanf(str, "%lf %d %d", &version, &format, &size) != 3) return 0;
- if(format){
- binary = true;
- Msg::Info("Mesh is in binary format");
- int one;
- if(fread(&one, sizeof(int), 1, fp) != 1) return 0;
- if(one != 1){
- swap = true;
- Msg::Info("Swapping bytes from binary file");
- }
- }
-
- }
- else if(!strncmp(&str[1], "PhysicalNames", 13)) {
-
- if(!fgets(str, sizeof(str), fp)) return 0;
- int numNames;
- if(sscanf(str, "%d", &numNames) != 1) return 0;
- for(int i = 0; i < numNames; i++) {
- int dim = -1, num;
- if(version > 2.0){
- if(fscanf(fp, "%d", &dim) != 1) return 0;
- }
- if(fscanf(fp, "%d", &num) != 1) return 0;
- if(!fgets(str, sizeof(str), fp)) return 0;
- std::string name = ExtractDoubleQuotedString(str, 256);
- if(name.size()) setPhysicalName(name, dim, num);
- }
-
- }
- else if(!strncmp(&str[1], "NO", 2) || !strncmp(&str[1], "Nodes", 5) ||
- !strncmp(&str[1], "ParametricNodes", 15)) {
-
- const bool parametric = !strncmp(&str[1], "ParametricNodes", 15);
- if(!fgets(str, sizeof(str), fp)) return 0;
- int numVertices;
- if(sscanf(str, "%d", &numVertices) != 1) return 0;
- Msg::Info("%d vertices", numVertices);
- Msg::ResetProgressMeter();
- vertexVector.clear();
- vertexMap.clear();
- int minVertex = numVertices + 1, maxVertex = -1;
- for(int i = 0; i < numVertices; i++) {
- int num;
- double xyz[3], uv[2];
- MVertex *newVertex = 0;
- if (!parametric){
- if(!binary){
- if (fscanf(fp, "%d %lf %lf %lf", &num, &xyz[0], &xyz[1], &xyz[2]) != 4)
- return 0;
- }
- else{
- if(fread(&num, sizeof(int), 1, fp) != 1) return 0;
- if(swap) SwapBytes((char*)&num, sizeof(int), 1);
- if(fread(xyz, sizeof(double), 3, fp) != 3) return 0;
- if(swap) SwapBytes((char*)xyz, sizeof(double), 3);
- }
- newVertex = new MVertex(xyz[0], xyz[1], xyz[2], 0, num);
- }
- else{
- int iClasDim, iClasTag;
- if(!binary){
- if (fscanf(fp, "%d %lf %lf %lf %d %d", &num, &xyz[0], &xyz[1], &xyz[2],
- &iClasDim, &iClasTag) != 6)
- return 0;
- }
- else{
- if(fread(&num, sizeof(int), 1, fp) != 1) return 0;
- if(swap) SwapBytes((char*)&num, sizeof(int), 1);
- if(fread(xyz, sizeof(double), 3, fp) != 3) return 0;
- if(swap) SwapBytes((char*)xyz, sizeof(double), 3);
- if(fread(&iClasDim, sizeof(int), 1, fp) != 1) return 0;
- if(swap) SwapBytes((char*)&iClasDim, sizeof(int), 1);
- if(fread(&iClasTag, sizeof(int), 1, fp) != 1) return 0;
- if(swap) SwapBytes((char*)&iClasTag, sizeof(int), 1);
- }
- if (iClasDim == 0){
- GVertex *gv = getVertexByTag(iClasTag);
- if (gv) gv->deleteMesh();
- newVertex = new MVertex(xyz[0], xyz[1], xyz[2], gv, num);
- }
- else if (iClasDim == 1){
- GEdge *ge = getEdgeByTag(iClasTag);
- if(!binary){
- if (fscanf(fp, "%lf", &uv[0]) != 1) return 0;
- }
- else{
- if(fread(uv, sizeof(double), 1, fp) != 1) return 0;
- if(swap) SwapBytes((char*)uv, sizeof(double), 1);
- }
- newVertex = new MEdgeVertex(xyz[0], xyz[1], xyz[2], ge, uv[0], -1.0, num);
- }
- else if (iClasDim == 2){
- GFace *gf = getFaceByTag(iClasTag);
- if(!binary){
- if (fscanf(fp, "%lf %lf", &uv[0], &uv[1]) != 2) return 0;
- }
- else{
- if(fread(uv, sizeof(double), 2, fp) != 2) return 0;
- if(swap) SwapBytes((char*)uv, sizeof(double), 2);
- }
- newVertex = new MFaceVertex(xyz[0], xyz[1], xyz[2], gf, uv[0], uv[1], num);
- }
- else if (iClasDim == 3){
- GRegion *gr = getRegionByTag(iClasTag);
- newVertex = new MVertex(xyz[0], xyz[1], xyz[2], gr, num);
- }
- }
- minVertex = std::min(minVertex, num);
- maxVertex = std::max(maxVertex, num);
- if(vertexMap.count(num))
- Msg::Warning("Skipping duplicate vertex %d", num);
- vertexMap[num] = newVertex;
- if(numVertices > 100000)
- Msg::ProgressMeter(i + 1, numVertices, "Reading nodes");
- }
- // If the vertex numbering is dense, tranfer the map into a
- // vector to speed up element creation
- if((int)vertexMap.size() == numVertices &&
- ((minVertex == 1 && maxVertex == numVertices) ||
- (minVertex == 0 && maxVertex == numVertices - 1))){
- Msg::Info("Vertex numbering is dense");
- vertexVector.resize(vertexMap.size() + 1);
- if(minVertex == 1)
- vertexVector[0] = 0;
- else
- vertexVector[numVertices] = 0;
- std::map<int, MVertex*>::const_iterator it = vertexMap.begin();
- for(; it != vertexMap.end(); ++it)
- vertexVector[it->first] = it->second;
- vertexMap.clear();
- }
-
- }
- else if(!strncmp(&str[1], "ELM", 3) || !strncmp(&str[1], "Elements", 8)) {
-
- if(!fgets(str, sizeof(str), fp)) return 0;
- int numElements;
- std::map<int, MElement*> parents;
- sscanf(str, "%d", &numElements);
- Msg::Info("%d elements", numElements);
- Msg::ResetProgressMeter();
- if(!binary){
- for(int i = 0; i < numElements; i++) {
- int num, type, physical = 0, elementary = 0, partition = 0, parentN = 0, numVertices;
- if(version <= 1.0){
- if(fscanf(fp, "%d %d %d %d %d", &num, &type, &physical, &elementary, &numVertices) != 5)
- return 0;
- if(numVertices != MElement::getInfoMSH(type)) return 0;
- }
- else{
- int numTags;
- if(fscanf(fp, "%d %d %d", &num, &type, &numTags) != 3) return 0;
- for(int j = 0; j < numTags; j++){
- int tag;
- if(fscanf(fp, "%d", &tag) != 1) return 0;
- if(j == 0) physical = tag;
- else if(j == 1) elementary = tag;
- else if(j == 2) partition = tag;
- else if(j == 3) parentN = tag;
- // ignore any other tags for now
- }
- if(!(numVertices = MElement::getInfoMSH(type))) {
- if(type != MSH_POLYG_ && type != MSH_POLYH_) return 0;
- if(fscanf(fp, "%d", &numVertices) != 1) return 0;
- }
- }
- int *indices = new int[numVertices];
- for(int j = 0; j < numVertices; j++)
- if(fscanf(fp, "%d", &indices[j]) != 1) return 0;
- std::vector<MVertex*> vertices;
- if(vertexVector.size()){
- if(!getVertices(numVertices, indices, vertexVector, vertices)) return 0;
- }
- else{
- if(!getVertices(numVertices, indices, vertexMap, vertices)) return 0;
- }
- MElement *e = createElementMSH(this, num, type, physical, elementary,
- partition, vertices, elements, physicals);
- // creation of interface elements (no interface element for a line)
- if(e->getNumPrimaryVertices() > 2) createInterfaceElementMSH(this, e, elementary, partition, vertices, elements, interfaces,physical);
- if(parentN) {
- MElement *parent = NULL;
- if(parents.find(parentN) == parents.end()){
- parent = getParent(parentN, e->getDim(), elements);
- if(parent) parents[parentN] = parent;
- else Msg::Error("Parent element %d not found!",parentN);
- }
- else parent = parents.find(parentN)->second;
- e->setParent(parent);
- }
- if(numElements > 100000)
- Msg::ProgressMeter(i + 1, numElements, "Reading elements");
- delete [] indices;
- }
- }
- else{
- int numElementsPartial = 0;
- while(numElementsPartial < numElements){
- int header[3];
- if(fread(header, sizeof(int), 3, fp) != 3) return 0;
- if(swap) SwapBytes((char*)header, sizeof(int), 3);
- int type = header[0];
- int numElms = header[1];
- int numTags = header[2];
- int numVertices = MElement::getInfoMSH(type);
- unsigned int n = 1 + numTags + numVertices;
- int *data = new int[n];
- for(int i = 0; i < numElms; i++) {
- if(fread(data, sizeof(int), n, fp) != n) return 0;
- if(swap) SwapBytes((char*)data, sizeof(int), n);
- int num = data[0];
- int physical = (numTags > 0) ? data[4 - numTags] : 0;
- int elementary = (numTags > 1) ? data[4 - numTags + 1] : 0;
- int partition = (numTags > 2) ? data[4 - numTags + 2] : 0;
- int *indices = &data[numTags + 1];
- std::vector<MVertex*> vertices;
- if(vertexVector.size()){
- if(!getVertices(numVertices, indices, vertexVector, vertices)) return 0;
- }
- else{
- if(!getVertices(numVertices, indices, vertexMap, vertices)) return 0;
- }
- createElementMSH(this, num, type, physical, elementary, partition,
- vertices, elements, physicals);
- if(numElements > 100000)
- Msg::ProgressMeter(numElementsPartial + i + 1, numElements,
- "Reading elements");
- }
- delete [] data;
- numElementsPartial += numElms;
- }
- }
- }
- else if(!strncmp(&str[1], "NodeData", 8)) {
-
- // there's some nodal post-processing data to read later on, so
- // cache the vertex indexing data
- if(vertexVector.size())
- _vertexVectorCache = vertexVector;
- else
- _vertexMapCache = vertexMap;
- postpro = true;
- break;
- }
- else if(!strncmp(&str[1], "ElementData", 11) ||
- !strncmp(&str[1], "ElementNodeData", 15)) {
-
- // there's some element post-processing data to read later on
- postpro = true;
- break;
- }
-
- do {
- if(!fgets(str, sizeof(str), fp) || feof(fp))
- break;
- } while(str[0] != '$');
- }
-
- // store the elements in their associated elementary entity. If the
- // entity does not exist, create a new (discrete) one.
- for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
- _storeElementsInEntities(elements[i]);
-
- // Store the interfaceElements when I understand this operation
- //_storeInterfaceElements(interelements);
-
- // associate the correct geometrical entity with each mesh vertex
- _associateEntityWithMeshVertices();
-
- // store the vertices in their associated geometrical entity
- if(vertexVector.size())
- _storeVerticesInEntities(vertexVector);
- else
- _storeVerticesInEntities(vertexMap);
-
- // store the physical tags
- for(int i = 0; i < 4; i++)
- _storePhysicalTagsInEntities(i, physicals[i]);
-
- fclose(fp);
-
- return postpro ? 2 : 1;
-}
-// Generate Interface Element on BoundaryCondition
-void GModelWithInterface::generateInterfaceElementsOnBoundary(const int &num_phys, std::vector<DGelasticField> elasticFields){
- //Boundary condition on the derivative of displacement
- std::vector<MVertex*> Mv;
- std::vector<MVertex*> vv;
- MElement *el=NULL;
- int ord,nedge,temp=0;
- bool flag_find;
- this->getMeshVerticesForPhysicalGroup(1,num_phys,Mv); // 1 because always 1D interface element for thin structures
- // Loop on nodes of Physical entities to find the associated element Optimize it because for each components of Mv loop on all elements (il faudrait les "noeuds interieurs du groupe physique)
- for(int j=0;j<Mv.size();j++){
- flag_find=false;
- for(int jj=0;jj<elasticFields.size();jj++){
- for (std::set<MElement*>::const_iterator it=elasticFields[jj].g->begin();it!=elasticFields[jj].g->end();++it){
- el = *it;
- // retrieve order and number of edge of element
- ord = el->getPolynomialOrder(); nedge = el->getNumEdges();
- // The interfaceElement is created only if Mv[j] is the first interior vertex (ie it is not a primary vertex) of the edge of the element
- // So only this vertex of element are tested
- for(int k=0;k<nedge;k++){
- temp = nedge + k*(ord-1);
- if (el->getVertex(temp) == Mv[j]){
- el->getEdgeVertices(k,vv);
- MInterfaceElement ie = MInterfaceElement(vv, 0, 0, el,el);
- boundinter.push_back(std::pair<int,MInterfaceElement>(jj,ie));
- vv.clear();
- flag_find=true;
- break;
- }
- else if(el->getVertex(k) == Mv[j]) {flag_find=true; break;} // Mv[j] is a primary vertices -->break // Remove it if Mv contains only the interior vertices
- }
- if(flag_find) break; // Element associated to the node is find -->break
- }
- if(flag_find) break; // Element associated to the node is find -->break
- }
- }
-}
diff --git a/contrib/cm3/GModelWithInterface.h b/contrib/cm3/GModelWithInterface.h
deleted file mode 100644
index ad07ff5171..0000000000
--- a/contrib/cm3/GModelWithInterface.h
+++ /dev/null
@@ -1,59 +0,0 @@
-//
-// C++ Interface: terms
-//
-// Description: Insertion of Interface Elements in GModel
-//
-//
-// Author: <Gauthier BECKER>, (C) 2010
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-
-
-#ifndef GMODELDG_H_
-#define GMODELDG_H_
-#include "GModel.h"
-#include "MInterfaceElement.h"
-#include "DgC0PlateSolver.h" // remove if elasticfields is not pass to generateInterfaceElementsOnBoundary
-#include "groupOfElements.h" // remove if elasticfields is not pass to generateInterfaceElementsOnBoundary
-class GModelWithInterface : public GModel{
-
- protected :
- // Add Interface element to the GModel I keep the MInterfaceElement because there is no need of a GInterfaceElement
- // I use vector because I don't know how to use set
- std::vector<std::pair<int,MInterfaceElement> > interfaces; // TODO make a pair to take into account more than 1 elasticfield
- std::vector<std::pair<int,MInterfaceElement> > boundinter; // A different vector is used for boundary element of the interface
-
- // Store interface elements when I understand this operation
- // void _storeInterfaceElements(std::vector<std::map<int,MInterfaceElement> > & ie );
-
- public :
- // build function
- GModelWithInterface() : GModel(){};
-
- // return interface
- std::vector<MInterfaceElement*> getInterface(const int num_field){
- std::vector<MInterfaceElement*> vie;
- for(int i=0;i<interfaces.size();i++)
- if(interfaces[i].first == num_field) vie.push_back(&interfaces[i].second);
- return vie;
- }
-
- // Function reading the input .msh file
- int readMSH(const std::string &name);
-
- // Function to generate interface element on boundary
- void generateInterfaceElementsOnBoundary(const int &num_phys, std::vector<DGelasticField> elasticFields);
-
- // Return the boundary interfaceElement linked to an elasticField
- std::vector<MInterfaceElement*> getBoundInterface(const int num_field){
- std::vector<MInterfaceElement*> vie;
- for(int i=0;i<boundinter.size();i++)
- if(boundinter[i].first == num_field) vie.push_back(&boundinter[i].second);
- return vie;
- }
-};
-
-
-#endif // GMODELDG_H_
diff --git a/contrib/cm3/LinearElasticShellHookeTensor.h b/contrib/cm3/LinearElasticShellHookeTensor.h
deleted file mode 100644
index e00a2dbbae..0000000000
--- a/contrib/cm3/LinearElasticShellHookeTensor.h
+++ /dev/null
@@ -1,75 +0,0 @@
-//
-// C++ Interface: terms
-//
-// Description: Hooke Tensor for shell (linear elastic)
-//
-//
-// Author: <Gauthier BECKER>, (C) 2010
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-
-# ifndef LINEARELASTICSHELLHOOKETENSOR_H_
-# define LINEARELASTICSHELLHOOKETENSOR_H_
-#include "Tensor4dim2.h"
-class LinearElasticShellHookeTensor : public Tensor4dim2{
-
- public :
- LinearElasticShellHookeTensor(){};
- ~LinearElasticShellHookeTensor(){};
- void set(const LocalBasis *lb,const double &C11,const double &nu){
- for(int i=0;i<2;i++)
- for(int j=0;j<2;j++)
- for(int ii=0;ii<2;ii++)
- for(int jj=0;jj<2;jj++)
- tensor[i][j][ii][jj] = C11*( nu*dot(lb->getphi0d(i),lb->getphi0d(j))*dot(lb->getphi0d(ii),lb->getphi0d(jj)) + 0.5*(1.-nu)*(dot(lb->getphi0d(i),lb->getphi0d(ii))*dot(lb->getphi0d(jj),lb->getphi0d(j)) + dot(lb->getphi0d(i),lb->getphi0d(jj))*dot(lb->getphi0d(ii),lb->getphi0d(j)) ) );
- }
-
- void set(const LocalBasis *lb,const double &nu){
- for(int i=0;i<2;i++)
- for(int j=0;j<2;j++)
- for(int ii=0;ii<2;ii++)
- for(int jj=0;jj<2;jj++)
- tensor[i][j][ii][jj] = ( nu*dot(lb->getphi0d(i),lb->getphi0d(j))*dot(lb->getphi0d(ii),lb->getphi0d(jj)) + 0.5*(1.-nu)*(dot(lb->getphi0d(i),lb->getphi0d(ii))*dot(lb->getphi0d(jj),lb->getphi0d(j)) + dot(lb->getphi0d(i),lb->getphi0d(jj))*dot(lb->getphi0d(ii),lb->getphi0d(j)) ) );
- }
-
- void hat(const LocalBasis *lb, const double C11, const double nu){
- double temp;
- for(int i=0;i<2;i++)
- for(int j=0;j<2;j++)
- for(int ii=0;ii<2;ii++)
- for(int jj=0;jj<2;jj++){
- temp =0.;
- for(int a=0;a<2;a++)
- for(int b=0;b<2;b++)
- for(int c=0;c<2;c++)
- for(int d=0;d<2;d++)
- temp += lb->getT(a,i)*lb->getT(b,j)*( nu*dot(lb->getphi0d(a),lb->getphi0d(b))*dot(lb->getphi0d(c),lb->getphi0d(d)) + 0.5*(1.-nu)*(dot(lb->getphi0d(a),lb->getphi0d(c))*dot(lb->getphi0d(d),lb->getphi0d(b)) + dot(lb->getphi0d(a),lb->getphi0d(d))*dot(lb->getphi0d(c),lb->getphi0d(b)) )) *lb->getT(c,ii)*lb->getT(d,jj);
- tensor[i][j][ii][jj]=C11*temp;
- }
- }
- void hat(const LocalBasis *lb,const LinearElasticShellHookeTensor *H){
- double temp;
- for(int i=0;i<2;i++)
- for(int j=0;j<2;j++)
- for(int ii=0;ii<2;ii++)
- for(int jj=0;jj<2;jj++){
- temp =0.;
- for(int a=0;a<2;a++)
- for(int b=0;b<2;b++)
- for(int c=0;c<2;c++)
- for(int d=0;d<2;d++)
- temp += lb->getT(a,i)*lb->getT(b,j)*H->get(a,b,c,d)*lb->getT(c,ii)*lb->getT(d,jj);
- this->tensor[i][j][ii][jj]=temp;
- }
- }
- void mean(const LinearElasticShellHookeTensor *Hp, const LinearElasticShellHookeTensor *Hm){
- for(int a=0;a<2;a++)
- for(int b=0;b<2;b++)
- for(int c=0;c<2;c++)
- for(int d=0;d<2;d++)
- tensor[a][b][c][d] = 0.5*(Hp->get(a,b,c,d)+Hm->get(a,b,c,d));
- }
-};
-#endif // LinearElasticShellHooketensor
diff --git a/contrib/cm3/LocalBasis.cpp b/contrib/cm3/LocalBasis.cpp
deleted file mode 100644
index e39f24e93a..0000000000
--- a/contrib/cm3/LocalBasis.cpp
+++ /dev/null
@@ -1,13 +0,0 @@
-//
-// C++ Interface: terms
-//
-// Description: Class of local basis for shells element
-//
-//
-// Author: <Gauthier BECKER>, (C) 2010
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-
-#include "LocalBasis.h"
diff --git a/contrib/cm3/LocalBasis.h b/contrib/cm3/LocalBasis.h
deleted file mode 100644
index 04ede3ff17..0000000000
--- a/contrib/cm3/LocalBasis.h
+++ /dev/null
@@ -1,174 +0,0 @@
-//
-// C++ Interface: terms
-//
-// Description: Define the localbasis of element for shells
-//
-//
-// Author: <Gauthier BECKER>, (C) 2010
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-
-
-#ifndef LOCALBASIS_H_
-#define LOCALBASIS_H_
-#include"SVector3.h"
-class LocalBasis{
- protected :
- std::vector<SVector3> phi0;
- std::vector<SVector3> phi0d;
- fullMatrix<double> T; // not used if interface element and bulk term --> create separate class ??
- fullMatrix<double> t; // not used if interface element and bulk term --> create separate class ??
- SVector3 t0;
- double detJ0;
-
- void setphiall(const double d){
- for(int i=0;i<2;i++)
- for(int j=0;j<3;j++){
- phi0[i](j)=d;
- }
- }
-
- public :
- LocalBasis(){
- phi0.reserve(2);
- phi0d.reserve(2);
- T.resize(2,2);
- t.resize(2,2);
- for(int i=0;i<2;i++){
- SVector3 temp(0.,0.,0.);
- phi0.push_back(temp);
- //phi0[i]=temp;
- phi0d.push_back(temp);
- }
- T.setAll(0.);
- t.setAll(0.);
- }
- ~LocalBasis(){}
-
- void set(MElement *ele, const std::vector<TensorialTraits<SVector3>::GradType> &Grads){
- const int nbFF = ele->getNumVertices();
- // Compute local basis vector
- this->setphiall(0.);
- for(int i=0;i<2;i++){
- for(int j=0;j<nbFF;j++){
- phi0[i](0) += Grads[j](0,i)*ele->getVertex(j)->x();
- phi0[i](1) += Grads[j](0,i)*ele->getVertex(j)->y();
- phi0[i](2) += Grads[j](0,i)*ele->getVertex(j)->z();
- }
- }
-
- // Compute normal and Jacobian
- t0 = crossprod(phi0[0],phi0[1]);
- detJ0=t0.norm();
- t0.normalize();
-
- // Compute dual basis vector
- double cos;
- for(int i=0;i<2;i++){
- phi0d[i] = crossprod(phi0[1-i],t0);
- cos = dot(phi0d[i],phi0[i]);
- phi0d[i] *= (1./cos);
- }
- }
-
- void set(MInterfaceElement *iele, const std::vector<TensorialTraits<SVector3>::GradType> &Grads, const SVector3 &t0p, const SVector3 &t0m){
- const int nbFF = iele->getNumVertices();
-
- // Compute of normal
- t0=t0p+t0m;
- assert(t0.norm()>1.e-8);
- t0.normalize();
-
- // Compute local basis vector
- this->setphiall(0.);
- for(int j=0;j<nbFF;j++){
- phi0[0](0) += Grads[j](0,0)*iele->getVertex(j)->x();
- phi0[0](1) += Grads[j](0,0)*iele->getVertex(j)->y();
- phi0[0](2) += Grads[j](0,0)*iele->getVertex(j)->z();
- }
- phi0[1] = crossprod(t0,phi0[0]);
- phi0[1].normalize();
-
- // Compute normal and Jacobian
- detJ0=crossprod(phi0[0],phi0[1]).norm();
-
- // Compute dual basis vector
- double cos;
- for(int i=0;i<2;i++){
- phi0d[i] = crossprod(phi0[1-i],t0);
- cos = dot(phi0d[i],phi0[i]);
- phi0d[i] *= (1./cos);
- }
- }
-
- void set(MInterfaceElement *iele, const std::vector<TensorialTraits<SVector3>::GradType> &Grads, const SVector3 &t0m){
- const int nbFF = iele->getNumVertices();
-
- // Compute of normal
- t0=t0m;
- assert(t0.norm()>1.e-8);
- t0.normalize();
-
- // Compute local basis vector
- this->setphiall(0.);
- for(int j=0;j<nbFF;j++){
- phi0[0](0) += Grads[j](0,0)*iele->getVertex(j)->x();
- phi0[0](1) += Grads[j](0,0)*iele->getVertex(j)->y();
- phi0[0](2) += Grads[j](0,0)*iele->getVertex(j)->z();
- }
- phi0[1] = crossprod(t0,phi0[0]);
- phi0[1].normalize();
-
- // Compute normal and Jacobian
- detJ0=crossprod(phi0[0],phi0[1]).norm();
-
- // Compute dual basis vector
- double cos;
- for(int i=0;i<2;i++){
- phi0d[i] = crossprod(phi0[1-i],t0);
- cos = dot(phi0d[i],phi0[i]);
- phi0d[i] *= (1./cos);
- }
- }
-
-
- // Push Forward Tensor (not used for interface element --> create a separate class ??
- void set_pushForward(const LocalBasis *lbs){
- // push forward tensor T and inverse push forward tensor t
- for(int i=0;i<2;i++)
- for(int j=0;j<2;j++){
- T(i,j)=dot(phi0[i],lbs->getphi0d(j));
- t(i,j)=dot(phi0d[j],lbs->getphi0(i));
- }
- }
- // get operation
- inline double getJacobian() const {return detJ0;}
- SVector3 gett0() const {return t0;}
- inline double gett0(const int i) const {return t0(i);}
- std::vector<SVector3> getphi0() const {return phi0;}
- std::vector<SVector3> getphi0d()const {return phi0d;}
- SVector3 getphi0(const int i) const {return phi0[i];}
- SVector3 getphi0d(const int i)const {return phi0d[i];}
- inline double getphi0(const int i,const int j) const {return phi0[i](j);}
- inline double getphi0d(const int i,const int j) const {return phi0d[i](j);}
- fullMatrix<double> getT() const {return T;}
- fullMatrix<double> gett() const {return t;}
- inline double getT(const int i, const int j) const {return T(i,j);}
- inline double gett(const int i, const int j) const {return t(i,j);}
-
- // Print
- void print(){
- printf("Basis phi0\n");
- printf("1 : %f %f %f\n",phi0[0](0),phi0[0](1),phi0[0](2));
- printf("2 : %f %f %f\n",phi0[1](0),phi0[1](1),phi0[1](2));
- printf("Normal : %f %f %f \n",t0(0),t0(1),t0(2));
- printf("Dual Basis phi0d\n");
- printf("1 : %f %f %f\n",phi0d[0](0),phi0d[0](1),phi0d[0](2));
- printf("2 : %f %f %f\n",phi0d[1](0),phi0d[1](1),phi0d[1](2));
- T.print("PushForward tensor");
- t.print("Inverse PushForward tensor");
- }
-};
-# endif // localbasis
diff --git a/contrib/cm3/MInterfaceElement.cpp b/contrib/cm3/MInterfaceElement.cpp
deleted file mode 100644
index 3c5e97520e..0000000000
--- a/contrib/cm3/MInterfaceElement.cpp
+++ /dev/null
@@ -1,13 +0,0 @@
-//
-// C++ Interface: terms
-//
-// Description: Class of interface element used for DG
-//
-//
-// Author: <Gauthier BECKER>, (C) 2010
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-
-#include "MInterfaceElement.h"
diff --git a/contrib/cm3/MInterfaceElement.h b/contrib/cm3/MInterfaceElement.h
deleted file mode 100644
index 6b30110e47..0000000000
--- a/contrib/cm3/MInterfaceElement.h
+++ /dev/null
@@ -1,167 +0,0 @@
-//
-// C++ Interface: terms
-//
-// Description: Class of interface element used for DG 2D only for the moment
-// thus the interface element is a line
-//
-// Author: <Gauthier BECKER>, (C) 2010
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-# ifndef _MINTERFACEELEMENT_H_
-# define _MINTERFACEELEMENT_H_
-
-#include "MLine.h"
-#include "SVector3.h"
-#include "MEdge.h"
-#include "quadratureRules.h"
-
-class MInterfaceElement : public MLineN{ // or don't derivate but in this case a vector with the vertices of interface element has to be save ??
- protected :
- // table of pointer on the two elements linked to the interface element
- MElement *_numElem[2];
- // edge's number linked to interface element of minus and plus element
- int _numEdge[2];
- // dir = true if the edge and the interface element are defined in the same sens and dir = false otherwise
- bool _dir[2];
-
- public :
- // Build function CG/DG case
- MInterfaceElement(MVertex *v0, MVertex *v1, std::vector<MVertex*> &v, int num = 0, int part = 0, MElement *e_minus = 0, MElement *e_plus = 0)
- : MLineN(v0, v1, v, num, part)
- {
- _numElem[0] = e_minus;
- _numElem[1] = e_plus;
-
- // Edge of element linked to interface element we identifie an interior point of the MLine (degree 2 min for plate) thus we used v[0]
- for(int jj=0;jj<2;jj++){
- int nopv = _numElem[jj]->getNumPrimaryVertices();
- std::vector<MVertex*> vv;
- for(int i = 0; i < nopv; i++){
- _numElem[jj]->getEdgeVertices(i,vv);
- for(unsigned int j = 2; j < vv.size(); j++){
- if(vv[j] == v[0]){
- _numEdge[jj] = i;
- if(v0 == vv[0]) _dir[jj] = true; // same orientation
- else _dir[jj] = false;
- }
- }
- }
- }
- }
- MInterfaceElement(std::vector<MVertex*> &v, int num = 0, int part = 0, MElement *e_minus = 0, MElement *e_plus = 0)
- : MLineN(v, num, part)
- {
- _numElem[0] = e_minus;
- _numElem[1] = e_plus;
-
- // Edge of element linked to interface element we identifie an interior point of the MLine (degree 2 min for plate) thus we used v[0]
- for(int jj=0;jj<2;jj++){
- int nopv = _numElem[jj]->getNumPrimaryVertices();
- std::vector<MVertex*> vv;
- for(int i = 0; i < nopv; i++){
- _numElem[jj]->getEdgeVertices(i,vv);
- for(unsigned int j = 2; j < vv.size(); j++){
- if(vv[j] == v[2]){ // v[2] because it is the first interior node
- _numEdge[jj] = i;
- if(v[0] == vv[0]) _dir[jj] = true; // same orientation
- else _dir[jj] = false;
- }
- }
- }
- }
- }
- // Destructor
- ~MInterfaceElement(){}
-
- // Give the number of the elements in a vector
- MElement ** getElem(){return _numElem;}
-
- // Give the number of minus 0 or plus 1 element
- MElement * getElem(int index){return _numElem[index];}
-
- // Get the u v value on element for a abscissa u on the interface element // TODO optmize by store in the class interface element the corresponding value (if many step must be compute once)??
- void getuvOnElem(const double u, double &uem, double &vem, double &uep, double &vep){ // w = 0 as no volume element are taken into account. The point is defined between u=-1 and u=1 on the interface element
- double ue=0.,ve=0.;
- for(int jj=0;jj<2;jj++){
- switch(_numElem[jj]->getType()){
- case TYPE_TRI :
- switch(_numEdge[jj]){
- case 0 :
- if(_dir[jj]) {ue = 0.5 * ( 1 + u ); ve = 0.;}
- else {ue = 0.5 * ( 1 - u ); ve = 0.;}
- break;
- case 1 :
- if(_dir[jj]) {ue = 0.5 * (1 - u) ; ve = 0.5 * ( 1 + u );}
- else {ue = 0.5 * (1 + u) ; ve = 0.5 * ( 1 - u );}
- break;
- case 2 :
- if(_dir[jj]) { ue = 0; ve = 0.5 * (1 - u);}
- else { ue = 0; ve = 0.5 * (1 + u);}
- break;
- }
- break;
- case TYPE_QUA :
- switch(_numEdge[jj]){
- case 0 :
- if(_dir[jj]) {ue = u; ve = -1.;}
- else {ue =-u; ve=-1;}
- break;
- case 1 :
- if(_dir[jj]) {ue =1.; ve = u;}
- else {ue = 1.; ve = -u;}
- break;
- case 2 :
- if(_dir[jj]) {ue = -u; ve = 1;}
- else {ue = u; ve = 1;}
- break;
- case 3 :
- if(_dir[jj]) {ue = -1; ve = -u;}
- else {ue = -1; ve = u;}
- break;
- }
- break;
- default : Msg::Error("The Method doesn't work for this type of element");
- }
- if(jj==0){uem=ue;vem=ve;}
- else {uep=ue;vep=ve;}
- }
- }
- // Compute the characteristic size of the side h_s = max_e (area_e/perimeter_e)
- double getCharacteristicSize(){
- //printf("minus\n");
- double cm = this->characSize(_numElem[0]);
- //printf("plus\n");
- double cp = this->characSize(_numElem[1]);
- return cm > cp ? cm : cp;
- }
-
- // This function can be defined as a method of MElement ??
- double characSize(MElement *e){
- // Compute the area of the element
- // jacobian value compute somewhere else --> Optimize it ?? (But change at each iteration)
- GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
- IntPt *GP;
- double perimeter = 0., Area = 0.;
- double jac[3][3];
- int npts=Integ_Bulk.getIntPoints(e,&GP);
- // Area = 0.;
- for( int i = 0; i < npts; i++){
- // Coordonate of Gauss' point i
- const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2];
- const double weight = GP[i].weight; const double detJ = e->getJacobian(u, v, w, jac); // Or compute jacobian with crossprod(phi0[0],phi0[1]) ??
- Area += weight * detJ;
- }
- // perimeter
- int nside = e->getNumEdges();
- for( int i = 0; i < nside; i++){
- // Distance between the two extremities
- MEdge edge = e->getEdge(i);
- perimeter += edge.getVertex(0)->distance(edge.getVertex(1));
- }
- return Area/perimeter;
- };
- bool getdir(const int i){return _dir[i];}
-};
-#endif // MInterfaceElement
diff --git a/contrib/cm3/Tensor4dim2.h b/contrib/cm3/Tensor4dim2.h
deleted file mode 100644
index af5627e118..0000000000
--- a/contrib/cm3/Tensor4dim2.h
+++ /dev/null
@@ -1,116 +0,0 @@
-//
-// C++ Interface: terms
-//
-// Description: 4th Tensor in 2D (16 components) uesd to compute Hooke tensor
-//
-//
-// Author: <Gauthier BECKER>, (C) 2010
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-
-# ifndef TENSOR4DIM2_H_
-# define TENSOR4DIM2_H_
-// 4th order tensor in 2D
-class Tensor4dim2{
- protected :
- double tensor[2][2][2][2];
- public :
- Tensor4dim2(){
- for(int a=0;a<2;a++)
- for(int b=0;b<2;b++)
- for(int c=0;c<2;c++)
- for(int d=0;d<2;d++)
- tensor[a][b][c][d]=0.;
- }
-
- Tensor4dim2(const Tensor4dim2 &_in){
- for(int i=0;i<2;i++)
- for(int j=0;j<2;j++)
- for(int ii=0;ii<2;ii++)
- for(int jj=0;jj<2;jj++){
- this->set(i,j,ii,jj,_in(i,j,ii,jj));
- }
- }
- ~Tensor4dim2(){}
-
- void set(const int a,const int b,const int c,const int d,const double dd){tensor[a][b][c][d]=dd;}
- void setAll(const double dd){
- for(int a=0;a<2;a++)
- for(int b=0;b<2;b++)
- for(int c=0;c<2;c++)
- for(int d=0;d<2;d++)
- tensor[a][b][c][d]=dd;
-
- }
-
- Tensor4dim2 & operator = (const Tensor4dim2 &_in){
- for(int i=0;i<2;i++)
- for(int j=0;j<2;j++)
- for(int ii=0;ii<2;ii++)
- for(int jj=0;jj<2;jj++){
- this->set(i,j,ii,jj,_in(i,j,ii,jj));
- }
- return *this;
- }
-
- double operator() (const int a,const int b,const int c,const int d) const {return tensor[a][b][c][d];}
- double get(const int a,const int b,const int c,const int d) const {return tensor[a][b][c][d];}
-
- Tensor4dim2 operator+ (Tensor4dim2 &H){
- Tensor4dim2 Hsum;
- double temp;
- for(int a=0;a<2;a++)
- for(int b=0;b<2;b++)
- for(int c=0;c<2;c++)
- for(int d=0;d<2;d++){
- temp=tensor[a][b][c][d]+H(a,b,c,d);
- Hsum.set(a,b,c,d,temp);
- }
- return Hsum;
- }
-
- Tensor4dim2 operator* (const double dd){
- Tensor4dim2 dH;
- double temp;
- for(int a=0;a<2;a++)
- for(int b=0;b<2;b++)
- for(int c=0;c<2;c++)
- for(int d=0;d<2;d++){
- temp=dd*tensor[a][b][c][d];
- dH.set(a,b,c,d,temp);
- }
- return dH;
- }
-
- void operator+= (const Tensor4dim2 &H){
- double temp;
- for(int a=0;a<2;a++)
- for(int b=0;b<2;b++)
- for(int c=0;c<2;c++)
- for(int d=0;d<2;d++){
- temp=tensor[a][b][c][d]+H(a,b,c,d);
- this->set(a,b,c,d,temp);
- }
- }
- void operator*= (const double dd){
- double temp;
- for(int a=0;a<2;a++)
- for(int b=0;b<2;b++)
- for(int c=0;c<2;c++)
- for(int d=0;d<2;d++){
- temp=dd*tensor[a][b][c][d];
- this->set(a,b,c,d,temp);
- }
- }
- void print(){
- for(int a=0;a<2;a++)
- for(int b=0;b<2;b++)
- for(int c=0;c<2;c++)
- for(int d=0;d<2;d++)
- printf("(%d,%d,%d,%d)= %f\n",a,b,c,d,tensor[a][b][c][d]);
-
- }
-};
-#endif // Tensor4dim2
diff --git a/contrib/cm3/mainDG.cpp b/contrib/cm3/mainDG.cpp
deleted file mode 100644
index 7f8f0382a0..0000000000
--- a/contrib/cm3/mainDG.cpp
+++ /dev/null
@@ -1,39 +0,0 @@
-#include "Gmsh.h"
-#include "DgC0PlateSolver.h"
-#include "PView.h"
-#include "PViewData.h"
-#include "../arc/highlevel.h"
-#include "groupOfElements.h"
-#include <iterator>
-#include <iostream>
-
-int main (int argc, char* argv[])
-{
-
- if (argc != 2){
- printf("usage : elasticity input_file_name\n");
- return -1;
- }
- GmshInitialize(argc, argv);
- // globals are still present in Gmsh
-
- // instanciate a solver
- DgC0PlateSolver mySolver (1000);
-
- // read some input file
- mySolver.readInputFile(argv[1]);
-
- // solve the problem
- mySolver.solve();
-
- PView *pv = mySolver.buildDisplacementView("displacement");
- pv->getData()->writeMSH("disp.msh", false);
- delete pv;
- //pv = mySolver.buildElasticEnergyView("elastic energy"); //error ??
- //pv->getData()->writeMSH("energ.msh", false);
- //delete pv;
-
- // stop gmsh
- GmshFinalize();
-
-}
--
GitLab