diff --git a/CMakeLists.txt b/CMakeLists.txt
index 62df3ea7a3cea927368bddfb315e3b89e7a2736b..1779b5ed4018e8c4829dc262c566c2a7b85973dd 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -18,6 +18,7 @@ project(gmsh CXX C)
option(ENABLE_ANN "Enable ANN to compute Approximate Nearest Neighbors" ON)
option(ENABLE_ARC "Enable ARC-related binary targets" OFF)
+option(ENABLE_CM3 "Enable CM3-related binary targets" OFF)
option(ENABLE_BLAS_LAPACK "Use BLAS/Lapack for basic linear algebra" ON)
option(ENABLE_CGNS "Enable CGNS mesh export" OFF)
option(ENABLE_CHACO "Enable Chaco mesh partitioner" ON)
@@ -79,7 +80,8 @@ set(GMSH_API
Mesh/meshGFaceDelaunayInsertion.h
Solver/dofManager.h Solver/femTerm.h Solver/laplaceTerm.h Solver/elasticityTerm.h
Solver/linearSystem.h Solver/linearSystemGMM.h Solver/linearSystemCSR.h
- Solver/linearSystemFull.h Solver/elasticitySolver.h
+ Solver/linearSystemFull.h Solver/elasticitySolver.h
+Solver/c0DgPlateSolver.h
Post/PView.h Post/PViewData.h Plugin/PluginManager.h
Graphics/drawContext.h
contrib/kbipack/gmp_normal_form.h contrib/kbipack/gmp_matrix.h
@@ -848,6 +850,9 @@ endif(CYGWIN)
if(ENABLE_ARC)
include(contrib/arc/CMakeLists.txt)
endif(ENABLE_ARC)
+if(ENABLE_CM3)
+ include(contrib/cm3/CMakeLists.txt)
+endif(ENABLE_CM3)
find_program(BISON bison)
find_program(FLEX flex)
diff --git a/Geo/GEntity.h b/Geo/GEntity.h
index 6dadf2b107094fc9762beb32d6914b5ba638b195..c42cbd36d126396103548cf1130d560973ac48af 100644
--- a/Geo/GEntity.h
+++ b/Geo/GEntity.h
@@ -99,7 +99,7 @@ class GEntity {
DiscreteSurface,
CompoundSurface,
Volume,
- DiscreteVolume,
+ DiscreteVolume,
CompoundVolume,
PartitionVertex,
PartitionCurve,
@@ -143,7 +143,7 @@ class GEntity {
"Discrete surface",
"Compound surface",
"Volume",
- "Discrete volume",
+ "Discrete volume",
"Compound Volume",
"Partition vertex",
"Partition curve",
@@ -220,7 +220,7 @@ class GEntity {
// get the oriented bounding box
virtual SOrientedBoundingBox getOBB() {return SOrientedBoundingBox(); }
-
+
// get/set the visibility flag
virtual char getVisibility();
virtual void setVisibility(char val, bool recursive=false){ _visible = val; }
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 1603cf92aa8c92d4d72108ff0e1f306baf28e132..45a117e0724404bbdbb40cd3635feaf78a531a43 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -33,7 +33,7 @@ class binding;
// A geometric model. The model is a "not yet" non-manifold B-Rep.
class GModel
{
- private:
+ protected:
// the name of the model
std::string _name;
@@ -53,7 +53,7 @@ class GModel
// ghost cell information (stores partitions for each element acting
// as a ghost cell)
std::multimap<MElement*, short> _ghostCells;
-
+
// an octree for fast mesh element lookup
Octree *_octree;
@@ -457,7 +457,7 @@ class GModel
static void registerBindings(binding *b);
// Store mesh elements of a chain in a new elementary and physical entity
- void storeChain(int dim, std::map<int, std::vector<MElement*> > &entityMap,
+ void storeChain(int dim, std::map<int, std::vector<MElement*> > &entityMap,
std::map<int, std::map<int, std::string> > &physicalMap)
{
_storeElementsInEntities(entityMap);
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 03745710437cd7d223c7bdf4d2d5bdf1c37b284e..b0c0db81d442b4de2caad31494b621a978241cc2 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -59,7 +59,7 @@ void GModel::_storePhysicalTagsInEntities(int dim,
if(ge){
std::map<int, std::string>::const_iterator it2 = it->second.begin();
for(; it2 != it->second.end(); ++it2){
- if(std::find(ge->physicals.begin(), ge->physicals.end(), it2->first) ==
+ if(std::find(ge->physicals.begin(), ge->physicals.end(), it2->first) ==
ge->physicals.end())
ge->physicals.push_back(it2->first);
}
@@ -67,7 +67,7 @@ void GModel::_storePhysicalTagsInEntities(int dim,
}
}
-static bool getVertices(int num, int *indices, std::map<int, MVertex*> &map,
+static bool getVertices(int num, int *indices, std::map<int, MVertex*> &map,
std::vector<MVertex*> &vertices)
{
for(int i = 0; i < num; i++){
@@ -131,11 +131,11 @@ static MElement *createElementMSH(GModel *m, int num, int typeMSH, int physical,
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;
}
@@ -195,14 +195,15 @@ int GModel::readMSH(const std::string &name)
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;
@@ -242,7 +243,7 @@ int GModel::readMSH(const std::string &name)
}
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;
@@ -253,13 +254,13 @@ int GModel::readMSH(const std::string &name)
vertexMap.clear();
int minVertex = numVertices + 1, maxVertex = -1;
for(int i = 0; i < numVertices; i++) {
- int num;
+ 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;
+ return 0;
}
else{
if(fread(&num, sizeof(int), 1, fp) != 1) return 0;
@@ -271,7 +272,7 @@ int GModel::readMSH(const std::string &name)
}
else{
int iClasDim, iClasTag;
- if(!binary){
+ if(!binary){
if (fscanf(fp, "%d %lf %lf %lf %d %d", &num, &xyz[0], &xyz[1], &xyz[2],
&iClasDim, &iClasTag) != 6)
return 0;
@@ -293,7 +294,7 @@ int GModel::readMSH(const std::string &name)
}
else if (iClasDim == 1){
GEdge *ge = getEdgeByTag(iClasTag);
- if(!binary){
+ if(!binary){
if (fscanf(fp, "%lf", &uv[0]) != 1) return 0;
}
else{
@@ -304,18 +305,18 @@ int GModel::readMSH(const std::string &name)
}
else if (iClasDim == 2){
GFace *gf = getFaceByTag(iClasTag);
- if(!binary){
+ 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);
+ newVertex = new MVertex(xyz[0], xyz[1], xyz[2], gr, num);
}
}
minVertex = std::min(minVertex, num);
@@ -323,17 +324,17 @@ int GModel::readMSH(const std::string &name)
if(vertexMap.count(num))
Msg::Warning("Skipping duplicate vertex %d", num);
vertexMap[num] = newVertex;
- if(numVertices > 100000)
+ if(numVertices > 100000)
Msg::ProgressMeter(i + 1, numVertices, "Reading nodes");
}
// If the vertex numbering is dense, transfer the map into a
// vector to speed up element creation
- if((int)vertexMap.size() == numVertices &&
+ 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)
+ if(minVertex == 1)
vertexVector[0] = 0;
else
vertexVector[numVertices] = 0;
@@ -408,7 +409,7 @@ int GModel::readMSH(const std::string &name)
else p = parents.find(parent)->second;
e->setParent(p);
}
- if(numElements > 100000)
+ if(numElements > 100000)
Msg::ProgressMeter(i + 1, numElements, "Reading elements");
delete [] indices;
}
@@ -448,7 +449,7 @@ int GModel::readMSH(const std::string &name)
for(int j = 0; j < numPartitions - 1; j++)
_ghostCells.insert(std::pair<MElement*, short>(e, -data[5 + j]));
if(numElements > 100000)
- Msg::ProgressMeter(numElementsPartial + i + 1, numElements,
+ Msg::ProgressMeter(numElementsPartial + i + 1, numElements,
"Reading elements");
}
delete [] data;
@@ -528,6 +529,7 @@ static void writeElementMSH(FILE *fp, GModel *model, T *ele, bool saveAll,
}
template<class T>
+
static void writeElementsMSH(FILE *fp, GModel *model, std::vector<T*> &ele,
bool saveAll, int saveSinglePartition, double version,
bool binary, int &num, int elementary,
@@ -635,7 +637,7 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
fprintf(fp, "$PhysicalNames\n");
fprintf(fp, "%d\n", numPhysicalNames());
for(piter it = firstPhysicalName(); it != lastPhysicalName(); it++)
- fprintf(fp, "%d %d \"%s\"\n", it->first.first, it->first.second,
+ fprintf(fp, "%d %d \"%s\"\n", it->first.first, it->first.second,
it->second.c_str());
fprintf(fp, "$EndPhysicalNames\n");
}
@@ -654,7 +656,7 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
getEntities(entities);
for(unsigned int i = 0; i < entities.size(); i++)
for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
- entities[i]->mesh_vertices[j]->writeMSH(fp, binary, saveParametric,
+ entities[i]->mesh_vertices[j]->writeMSH(fp, binary, saveParametric,
scalingFactor);
if(binary) fprintf(fp, "\n");
@@ -782,7 +784,7 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
writeElementsMSH(fp, this, (*it)->polyhedra, saveAll, saveSinglePartition,
version, binary, num, (*it)->tag(), (*it)->physicals,
&parentsNum);
-
+
if(binary) fprintf(fp, "\n");
if(version >= 2.0){
@@ -1148,7 +1150,7 @@ int GModel::writeDistanceMSH(const std::string &name, double version, bool binar
int GModel::writePOS(const std::string &name, bool printElementary,
bool printElementNumber, bool printGamma, bool printEta,
- bool printRho, bool printDisto, bool saveAll,
+ bool printRho, bool printDisto, bool saveAll,
double scalingFactor)
{
FILE *fp = fopen(name.c_str(), "w");
@@ -1383,7 +1385,7 @@ int GModel::writeSTL(const std::string &name, bool binary, bool saveAll,
}
fwrite(&nfacets, sizeof(unsigned int), 1, fp);
}
-
+
for(fiter it = firstFace(); it != lastFace(); ++it) {
if(saveAll || (*it)->physicals.size()){
for(unsigned int i = 0; i < (*it)->triangles.size(); i++)
@@ -1395,7 +1397,7 @@ int GModel::writeSTL(const std::string &name, bool binary, bool saveAll,
if(!binary)
fprintf(fp, "endsolid Created by Gmsh\n");
-
+
fclose(fp);
return 1;
}
@@ -1406,9 +1408,9 @@ static int skipUntil(FILE *fp, const char *key)
strcpy(key_bracket, key);
strcat(key_bracket, "[");
while(fscanf(fp, "%s", str)){
- if(!strcmp(str, key)){
+ if(!strcmp(str, key)){
while(!feof(fp) && fgetc(fp) != '['){}
- return 1;
+ return 1;
}
if(!strcmp(str, key_bracket)){
return 1;
@@ -1432,7 +1434,7 @@ static int readVerticesVRML(FILE *fp, std::vector<MVertex*> &vertexVector,
}
static int readElementsVRML(FILE *fp, std::vector<MVertex*> &vertexVector, int region,
- std::map<int, std::vector<MElement*> > elements[3],
+ std::map<int, std::vector<MElement*> > elements[3],
bool strips=false)
{
int i;
@@ -1494,7 +1496,7 @@ static int readElementsVRML(FILE *fp, std::vector<MVertex*> &vertexVector, int r
Msg::Error("Prematured end of VRML file");
return 0;
}
- Msg::Info("%d elements", elements[0][region].size() +
+ Msg::Info("%d elements", elements[0][region].size() +
elements[1][region].size() + elements[2][region].size());
return 1;
}
@@ -1560,7 +1562,7 @@ int GModel::readVRML(const std::string &name)
}
}
- for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
+ for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
_storeElementsInEntities(elements[i]);
_associateEntityWithMeshVertices();
_storeVerticesInEntities(allVertexVector);
@@ -1587,13 +1589,13 @@ int GModel::writeVRML(const std::string &name, bool saveAll, double scalingFacto
fprintf(fp, " point [\n");
for(viter it = firstVertex(); it != lastVertex(); ++it)
- for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
+ for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
(*it)->mesh_vertices[i]->writeVRML(fp, scalingFactor);
for(eiter it = firstEdge(); it != lastEdge(); ++it)
for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
(*it)->mesh_vertices[i]->writeVRML(fp, scalingFactor);
for(fiter it = firstFace(); it != lastFace(); ++it)
- for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
+ for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
(*it)->mesh_vertices[i]->writeVRML(fp, scalingFactor);
fprintf(fp, " ]\n");
@@ -1666,7 +1668,7 @@ int GModel::readUNV(const std::string &name)
if(strlen(buffer) < 3) continue; // possible line ending after last fscanf
if(!strncmp(buffer, " -1", 6)) break;
int num, type, elementary, physical, color, numNodes;
- if(sscanf(buffer, "%d %d %d %d %d %d", &num, &type, &elementary, &physical,
+ if(sscanf(buffer, "%d %d %d %d %d %d", &num, &type, &elementary, &physical,
&color, &numNodes) != 6) break;
if(elementary < 0) elementary = 1;
if(physical < 0) physical = 0;
@@ -1695,17 +1697,17 @@ int GModel::readUNV(const std::string &name)
(new MLine3(vertices[0], vertices[2], vertices[1], num));
dim = 1;
break;
- case 41: case 51: case 61: case 74: case 81: case 91:
+ case 41: case 51: case 61: case 74: case 81: case 91:
elements[1][elementary].push_back(new MTriangle(vertices, num));
dim = 2;
break;
- case 42: case 52: case 62: case 72: case 82: case 92:
+ case 42: case 52: case 62: case 72: case 82: case 92:
elements[1][elementary].push_back
- (new MTriangle6(vertices[0], vertices[2], vertices[4], vertices[1],
+ (new MTriangle6(vertices[0], vertices[2], vertices[4], vertices[1],
vertices[3], vertices[5], num));
dim = 2;
break;
- case 44: case 54: case 64: case 71: case 84: case 94:
+ case 44: case 54: case 64: case 71: case 84: case 94:
elements[2][elementary].push_back(new MQuadrangle(vertices, num));
dim = 2;
break;
@@ -1718,64 +1720,64 @@ int GModel::readUNV(const std::string &name)
break;
case 111:
elements[3][elementary].push_back(new MTetrahedron(vertices, num));
- dim = 3;
+ dim = 3;
break;
- case 118:
+ case 118:
elements[3][elementary].push_back
(new MTetrahedron10(vertices[0], vertices[2], vertices[4], vertices[9],
vertices[1], vertices[3], vertices[5], vertices[6],
vertices[8], vertices[7], num));
dim = 3;
break;
- case 104: case 115:
+ case 104: case 115:
elements[4][elementary].push_back(new MHexahedron(vertices, num));
- dim = 3;
+ dim = 3;
break;
case 105: case 116:
elements[4][elementary].push_back
- (new MHexahedron20(vertices[0], vertices[2], vertices[4], vertices[6],
+ (new MHexahedron20(vertices[0], vertices[2], vertices[4], vertices[6],
vertices[12], vertices[14], vertices[16], vertices[18],
vertices[1], vertices[7], vertices[8], vertices[3],
vertices[9], vertices[5], vertices[10], vertices[11],
- vertices[13], vertices[19], vertices[15], vertices[17],
+ vertices[13], vertices[19], vertices[15], vertices[17],
num));
- dim = 3;
+ dim = 3;
break;
case 101: case 112:
elements[5][elementary].push_back(new MPrism(vertices, num));
- dim = 3;
+ dim = 3;
break;
case 102: case 113:
elements[5][elementary].push_back
- (new MPrism15(vertices[0], vertices[2], vertices[4], vertices[9],
- vertices[11], vertices[13], vertices[1], vertices[5],
+ (new MPrism15(vertices[0], vertices[2], vertices[4], vertices[9],
+ vertices[11], vertices[13], vertices[1], vertices[5],
vertices[6], vertices[3], vertices[7], vertices[8],
vertices[10], vertices[14], vertices[12], num));
dim = 3;
break;
}
-
- if(dim >= 0 && physical && (!physicals[dim].count(elementary) ||
+
+ if(dim >= 0 && physical && (!physicals[dim].count(elementary) ||
!physicals[dim][elementary].count(physical)))
physicals[dim][elementary][physical] = "unnamed";
}
}
}
}
-
- for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
+
+ for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
_storeElementsInEntities(elements[i]);
_associateEntityWithMeshVertices();
_storeVerticesInEntities(vertexMap);
- for(int i = 0; i < 4; i++)
+ for(int i = 0; i < 4; i++)
_storePhysicalTagsInEntities(i, physicals[i]);
fclose(fp);
return 1;
}
-int GModel::writeUNV(const std::string &name, bool saveAll, bool saveGroupsOfNodes,
+int GModel::writeUNV(const std::string &name, bool saveAll, bool saveGroupsOfNodes,
double scalingFactor)
{
FILE *fp = fopen(name.c_str(), "w");
@@ -1795,7 +1797,7 @@ int GModel::writeUNV(const std::string &name, bool saveAll, bool saveGroupsOfNod
fprintf(fp, "%6d\n", -1);
fprintf(fp, "%6d\n", 2411);
for(unsigned int i = 0; i < entities.size(); i++)
- for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
+ for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
entities[i]->mesh_vertices[j]->writeUNV(fp, scalingFactor);
fprintf(fp, "%6d\n", -1);
@@ -1834,7 +1836,7 @@ int GModel::writeUNV(const std::string &name, bool saveAll, bool saveGroupsOfNod
nodes.insert(e->getVertex(k));
}
}
- fprintf(fp, "%10d%10d%10d%10d%10d%10d%10d%10d\n",
+ fprintf(fp, "%10d%10d%10d%10d%10d%10d%10d%10d\n",
gr, 0, 0, 0, 0, 0, 0, (int)nodes.size());
fprintf(fp, "PERMANENT GROUP%d\n", gr++);
int row = 0;
@@ -1963,7 +1965,7 @@ int GModel::readMESH(const std::string &name)
}
}
- for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
+ for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
_storeElementsInEntities(elements[i]);
_associateEntityWithMeshVertices();
_storeVerticesInEntities(vertexVector);
@@ -1972,7 +1974,7 @@ int GModel::readMESH(const std::string &name)
return 1;
}
-int GModel::writeMESH(const std::string &name, int elementTagType,
+int GModel::writeMESH(const std::string &name, int elementTagType,
bool saveAll, double scalingFactor)
{
FILE *fp = fopen(name.c_str(), "w");
@@ -1994,9 +1996,9 @@ int GModel::writeMESH(const std::string &name, int elementTagType,
std::vector<GEntity*> entities;
getEntities(entities);
for(unsigned int i = 0; i < entities.size(); i++)
- for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
+ for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
entities[i]->mesh_vertices[j]->writeMESH(fp, scalingFactor);
-
+
int numEdges = 0, numTriangles = 0, numQuadrangles = 0, numTetrahedra = 0;
for(eiter it = firstEdge(); it != lastEdge(); ++it){
if(saveAll || (*it)->physicals.size()){
@@ -2022,7 +2024,7 @@ int GModel::writeMESH(const std::string &name, int elementTagType,
int numPhys = (*it)->physicals.size();
if(saveAll || numPhys){
for(unsigned int i = 0; i < (*it)->lines.size(); i++)
- (*it)->lines[i]->writeMESH(fp, elementTagType, (*it)->tag(),
+ (*it)->lines[i]->writeMESH(fp, elementTagType, (*it)->tag(),
numPhys ? (*it)->physicals[0] : 0);
}
}
@@ -2034,7 +2036,7 @@ int GModel::writeMESH(const std::string &name, int elementTagType,
int numPhys = (*it)->physicals.size();
if(saveAll || numPhys){
for(unsigned int i = 0; i < (*it)->triangles.size(); i++)
- (*it)->triangles[i]->writeMESH(fp, elementTagType, (*it)->tag(),
+ (*it)->triangles[i]->writeMESH(fp, elementTagType, (*it)->tag(),
numPhys ? (*it)->physicals[0] : 0);
}
}
@@ -2065,12 +2067,12 @@ int GModel::writeMESH(const std::string &name, int elementTagType,
}
fprintf(fp, " End\n");
-
+
fclose(fp);
return 1;
}
-int GModel::writeFEA(const std::string &name, int elementTagType,
+int GModel::writeFEA(const std::string &name, int elementTagType,
bool saveAll, double scalingFactor)
{
FILE *fp = fopen(name.c_str(), "w");
@@ -2083,7 +2085,7 @@ int GModel::writeFEA(const std::string &name, int elementTagType,
int numVertices = indexMeshVertices(saveAll), num2D = 0, num3D = 0;
for(fiter it = firstFace(); it != lastFace(); ++it)
- if(saveAll || (*it)->physicals.size())
+ if(saveAll || (*it)->physicals.size())
num2D += (*it)->getNumMeshElements();
for(riter it = firstRegion(); it != lastRegion(); ++it)
if(saveAll || (*it)->physicals.size())
@@ -2094,12 +2096,12 @@ int GModel::writeFEA(const std::string &name, int elementTagType,
fprintf(fp,"%d %d %d\n", numVertices, num2D, num3D);
else
fprintf(fp,"%d %d\n", numVertices, num2D ? num2D : num3D);
-
+
std::vector<GEntity*> entities;
getEntities(entities);
for(unsigned int i = 0; i < entities.size(); i++)
- for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
+ for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
if(entities[i]->mesh_vertices[j]->getIndex() >= 0)
fprintf(fp,"%d %g %g %g\n", entities[i]->mesh_vertices[j]->getIndex(),
entities[i]->mesh_vertices[j]->x() * scalingFactor,
@@ -2112,7 +2114,7 @@ int GModel::writeFEA(const std::string &name, int elementTagType,
if(saveAll || numPhys)
for(unsigned int i = 0; i < (*it)->getNumMeshElements(); i++)
(*it)->getMeshElement(i)->writeFEA
- (fp, elementTagType, iElement++, (*it)->tag(),
+ (fp, elementTagType, iElement++, (*it)->tag(),
numPhys ? (*it)->physicals[0] : 0);
}
@@ -2121,10 +2123,10 @@ int GModel::writeFEA(const std::string &name, int elementTagType,
if(saveAll || numPhys)
for(unsigned int i = 0; i < (*it)->getNumMeshElements(); i++)
(*it)->getMeshElement(i)->writeFEA
- (fp, elementTagType, iElement++, (*it)->tag(),
+ (fp, elementTagType, iElement++, (*it)->tag(),
numPhys ? (*it)->physicals[0] : 0);
}
-
+
fclose(fp);
return 1;
}
@@ -2145,8 +2147,8 @@ static double atofBDF(char *str)
if(str[i] == 'E' || str[i] == 'e') {
return atof(str);
}
- else if(str[i] == 'D' || str[i] == 'd'){
- str[i] = 'E';
+ else if(str[i] == 'D' || str[i] == 'd'){
+ str[i] = 'E';
return atof(str);
}
}
@@ -2164,7 +2166,7 @@ static double atofBDF(char *str)
return atof(tmp);
}
-static int readVertexBDF(FILE *fp, char *buffer, int keySize,
+static int readVertexBDF(FILE *fp, char *buffer, int keySize,
int *num, double *x, double *y, double *z)
{
char tmp[5][32];
@@ -2183,9 +2185,9 @@ static int readVertexBDF(FILE *fp, char *buffer, int keySize,
case 1: // small field
for(int i = 0; i < 5; i++) tmp[i][8] = '\0';
strncpy(tmp[0], &buffer[8], 8);
- strncpy(tmp[2], &buffer[24], 8);
- strncpy(tmp[3], &buffer[32], 8);
- strncpy(tmp[4], &buffer[40], 8);
+ strncpy(tmp[2], &buffer[24], 8);
+ strncpy(tmp[3], &buffer[32], 8);
+ strncpy(tmp[4], &buffer[40], 8);
break;
case 2: // long field
for(int i = 0; i < 5; i++) tmp[i][16] = '\0';
@@ -2232,7 +2234,7 @@ static void readLineBDF(char *buffer, int format, std::vector<char*> &fields)
}
}
-static int readElementBDF(FILE *fp, char *buffer, int keySize, int numVertices,
+static int readElementBDF(FILE *fp, char *buffer, int keySize, int numVertices,
int &num, int ®ion, std::vector<MVertex*> &vertices,
std::map<int, MVertex*> &vertexMap)
{
@@ -2244,14 +2246,14 @@ static int readElementBDF(FILE *fp, char *buffer, int keySize, int numVertices,
readLineBDF(buffer, format, fields);
- if(((int)fields.size() - 2 < abs(numVertices)) ||
+ if(((int)fields.size() - 2 < abs(numVertices)) ||
(numVertices < 0 && (fields.size() == 9))){
if(fields.size() == 9) fields.pop_back();
if(!fgets(buffer2, sizeof(buffer2), fp)) return 0;
readLineBDF(buffer2, format, fields);
}
- if(((int)fields.size() - 2 < abs(numVertices)) ||
+ if(((int)fields.size() - 2 < abs(numVertices)) ||
(numVertices < 0 && (fields.size() == 17))){
if(fields.size() == 17) fields.pop_back();
if(!fgets(buffer3, sizeof(buffer3), fp)) return 0;
@@ -2356,8 +2358,8 @@ int GModel::readBDF(const std::string &name)
if(readElementBDF(fp, buffer, 6, -4, num, region, vertices, vertexMap)){
if(vertices.size() == 10)
elements[3][region].push_back
- (new MTetrahedron10(vertices[0], vertices[1], vertices[2], vertices[3],
- vertices[4], vertices[5], vertices[6], vertices[7],
+ (new MTetrahedron10(vertices[0], vertices[1], vertices[2], vertices[3],
+ vertices[4], vertices[5], vertices[6], vertices[7],
vertices[9], vertices[8], num));
else
elements[3][region].push_back(new MTetrahedron(vertices, num));
@@ -2367,11 +2369,11 @@ int GModel::readBDF(const std::string &name)
if(readElementBDF(fp, buffer, 5, -8, num, region, vertices, vertexMap)){
if(vertices.size() == 20)
elements[4][region].push_back
- (new MHexahedron20(vertices[0], vertices[1], vertices[2], vertices[3],
- vertices[4], vertices[5], vertices[6], vertices[7],
- vertices[8], vertices[11], vertices[12], vertices[9],
- vertices[13], vertices[10], vertices[14], vertices[15],
- vertices[16], vertices[19], vertices[17], vertices[18],
+ (new MHexahedron20(vertices[0], vertices[1], vertices[2], vertices[3],
+ vertices[4], vertices[5], vertices[6], vertices[7],
+ vertices[8], vertices[11], vertices[12], vertices[9],
+ vertices[13], vertices[10], vertices[14], vertices[15],
+ vertices[16], vertices[19], vertices[17], vertices[18],
num));
else
elements[4][region].push_back(new MHexahedron(vertices, num));
@@ -2395,8 +2397,8 @@ int GModel::readBDF(const std::string &name)
}
}
}
-
- for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
+
+ for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
_storeElementsInEntities(elements[i]);
_associateEntityWithMeshVertices();
_storeVerticesInEntities(vertexMap);
@@ -2405,7 +2407,7 @@ int GModel::readBDF(const std::string &name)
return 1;
}
-int GModel::writeBDF(const std::string &name, int format, int elementTagType,
+int GModel::writeBDF(const std::string &name, int format, int elementTagType,
bool saveAll, double scalingFactor)
{
FILE *fp = fopen(name.c_str(), "w");
@@ -2425,7 +2427,7 @@ int GModel::writeBDF(const std::string &name, int format, int elementTagType,
// nodes
for(unsigned int i = 0; i < entities.size(); i++)
- for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
+ for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
entities[i]->mesh_vertices[j]->writeBDF(fp, format, scalingFactor);
// elements
@@ -2434,12 +2436,12 @@ int GModel::writeBDF(const std::string &name, int format, int elementTagType,
int numPhys = entities[i]->physicals.size();
if(saveAll || numPhys)
entities[i]->getMeshElement(j)->writeBDF
- (fp, format, elementTagType, entities[i]->tag(),
+ (fp, format, elementTagType, entities[i]->tag(),
numPhys ? entities[i]->physicals[0] : 0);
}
fprintf(fp, "ENDDATA\n");
-
+
fclose(fp);
return 1;
}
@@ -2537,8 +2539,8 @@ int GModel::readP3D(const std::string &name)
gr->transfinite_vertices[i + 1][j + 1][k + 1],
gr->transfinite_vertices[i ][j + 1][k + 1]));
}
- }
-
+ }
+
fclose(fp);
return 1;
}
@@ -2555,17 +2557,17 @@ int GModel::writeP3D(const std::string &name, bool saveAll, double scalingFactor
std::vector<GFace*> faces;
for(fiter it = firstFace(); it != lastFace(); ++it)
- if((*it)->transfinite_vertices.size() &&
+ if((*it)->transfinite_vertices.size() &&
(*it)->transfinite_vertices[0].size() &&
((*it)->physicals.size() || saveAll)) faces.push_back(*it);
std::vector<GRegion*> regions;
for(riter it = firstRegion(); it != lastRegion(); ++it)
- if((*it)->transfinite_vertices.size() &&
- (*it)->transfinite_vertices[0].size() &&
- (*it)->transfinite_vertices[0][0].size() &&
+ if((*it)->transfinite_vertices.size() &&
+ (*it)->transfinite_vertices[0].size() &&
+ (*it)->transfinite_vertices[0][0].size() &&
((*it)->physicals.size() || saveAll)) regions.push_back(*it);
-
+
if(faces.empty() && regions.empty()){
Msg::Warning("No structured grids to save");
fclose(fp);
@@ -2573,18 +2575,18 @@ int GModel::writeP3D(const std::string &name, bool saveAll, double scalingFactor
}
fprintf(fp, "%d\n", (int)(faces.size() + regions.size()));
-
+
for(unsigned int i = 0; i < faces.size(); i++)
- fprintf(fp, "%d %d 1\n",
+ fprintf(fp, "%d %d 1\n",
(int)faces[i]->transfinite_vertices.size(),
(int)faces[i]->transfinite_vertices[0].size());
for(unsigned int i = 0; i < regions.size(); i++)
- fprintf(fp, "%d %d %d\n",
+ fprintf(fp, "%d %d %d\n",
(int)regions[i]->transfinite_vertices.size(),
(int)regions[i]->transfinite_vertices[0].size(),
(int)regions[i]->transfinite_vertices[0][0].size());
-
+
for(unsigned int i = 0; i < faces.size(); i++){
GFace *gf = faces[i];
for(int coord = 0; coord < 3; coord++){
@@ -2598,7 +2600,7 @@ int GModel::writeP3D(const std::string &name, bool saveAll, double scalingFactor
}
}
}
-
+
for(unsigned int i = 0; i < regions.size(); i++){
GRegion *gr = regions[i];
for(int coord = 0; coord < 3; coord++){
@@ -2614,7 +2616,7 @@ int GModel::writeP3D(const std::string &name, bool saveAll, double scalingFactor
}
}
}
-
+
fclose(fp);
return 1;
}
@@ -2649,10 +2651,10 @@ int GModel::writeVTK(const std::string &name, bool binary, bool saveAll,
// write mesh vertices
fprintf(fp, "POINTS %d double\n", numVertices);
for(unsigned int i = 0; i < entities.size(); i++)
- for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
+ for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
entities[i]->mesh_vertices[j]->writeVTK(fp, binary, scalingFactor, bigEndian);
fprintf(fp, "\n");
-
+
// loop over all elements we need to save and count vertices
int numElements = 0, totalNumInt = 0;
for(unsigned int i = 0; i < entities.size(); i++){
@@ -2665,7 +2667,7 @@ int GModel::writeVTK(const std::string &name, bool binary, bool saveAll,
}
}
}
-
+
// print vertex indices in ascii or binary
fprintf(fp, "CELLS %d %d\n", numElements, totalNumInt);
for(unsigned int i = 0; i < entities.size(); i++){
@@ -2677,8 +2679,8 @@ int GModel::writeVTK(const std::string &name, bool binary, bool saveAll,
}
}
fprintf(fp, "\n");
-
- // print element types in ascii or binary
+
+ // print element types in ascii or binary
fprintf(fp, "CELL_TYPES %d\n", numElements);
for(unsigned int i = 0; i < entities.size(); i++){
if(entities[i]->physicals.size() || saveAll){
@@ -2696,7 +2698,7 @@ int GModel::writeVTK(const std::string &name, bool binary, bool saveAll,
}
}
}
-
+
fclose(fp);
return 1;
}
@@ -2708,23 +2710,23 @@ int GModel::readVTK(const std::string &name, bool bigEndian)
Msg::Error("Unable to open file '%s'", name.c_str());
return 0;
}
-
+
char buffer[256], buffer2[256];
-
+
if(!fgets(buffer, sizeof(buffer), fp)) return 0; // version line
if(!fgets(buffer, sizeof(buffer), fp)) return 0; // title
-
+
if(fscanf(fp, "%s", buffer) != 1) // ASCII or BINARY
Msg::Error("Failed reading buffer");
bool binary = false;
if(!strcmp(buffer, "BINARY")) binary = true;
-
+
if(fscanf(fp, "%s %s", buffer, buffer2) != 2) return 0;
if(strcmp(buffer, "DATASET") || strcmp(buffer2, "UNSTRUCTURED_GRID")){
Msg::Error("VTK reader can only read unstructured datasets");
return 0;
}
-
+
// read mesh vertices
int numVertices;
if(fscanf(fp, "%s %d %s\n", buffer, &numVertices, buffer2) != 3) return 0;
@@ -2818,17 +2820,17 @@ int GModel::readVTK(const std::string &name, bool bigEndian)
case 12: elements[5][1].push_back(new MHexahedron(cells[i])); break;
case 13: elements[6][1].push_back(new MPrism(cells[i])); break;
case 14: elements[7][1].push_back(new MPyramid(cells[i])); break;
- default:
+ default:
Msg::Error("Unknown type of cell %d", type);
break;
}
- }
-
- for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
+ }
+
+ for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
_storeElementsInEntities(elements[i]);
_associateEntityWithMeshVertices();
_storeVerticesInEntities(vertices);
-
+
fclose(fp);
return 1;
}
@@ -2846,17 +2848,17 @@ int GModel::readDIFF(const std::string &name)
std::map<int, std::map<int, std::string> > physicals[4];
std::map<int, MVertex*> vertexMap;
std::vector<MVertex*> vertexVector;
-
+
while(1) {
while(strstr(str, "Number of space dim. =") == NULL){
if(!fgets(str, sizeof(str), fp) || feof(fp))
break;
}
-
+
int dim;
if(sscanf(str, "%*s %*s %*s %*s %*s %d", &dim) != 1) return 0;
- Msg::Info("dimension %d", dim);
+ Msg::Info("dimension %d", dim);
int numElements;
if(!fgets(str, sizeof(str), fp) || feof(fp)) return 0;
@@ -2865,7 +2867,7 @@ int GModel::readDIFF(const std::string &name)
break;
}
if(sscanf(str, "%*s %*s %*s %*s %d", &numElements) != 1) return 0;
- Msg::Info("%d elements", numElements);
+ Msg::Info("%d elements", numElements);
int numVertices;
if(!fgets(str, sizeof(str), fp) || feof(fp)) return 0;
@@ -2874,7 +2876,7 @@ int GModel::readDIFF(const std::string &name)
break;
}
if(sscanf(str, "%*s %*s %*s %*s %d", &numVertices) != 1) return 0;
- Msg::Info("%d vertices", numVertices);
+ Msg::Info("%d vertices", numVertices);
int numVerticesPerElement;
if(!fgets(str, sizeof(str), fp)||feof(fp)) return 0;
@@ -2884,7 +2886,7 @@ int GModel::readDIFF(const std::string &name)
}
if(sscanf(str, "%*s %*s %*s %*s %*s %*s %*s %d", &numVerticesPerElement) != 1)
return 0;
- Msg::Info("numVerticesPerElement %d",numVerticesPerElement);
+ Msg::Info("numVerticesPerElement %d",numVerticesPerElement);
bool several_subdomains;
if(!fgets(str, sizeof(str), fp) || feof(fp)) return 0;
@@ -2896,8 +2898,8 @@ int GModel::readDIFF(const std::string &name)
several_subdomains = true;
else
several_subdomains = false;
- Msg::Info("several_subdomains %x", several_subdomains);
-
+ Msg::Info("several_subdomains %x", several_subdomains);
+
int nbi;
std::vector<int> bi;
if(!fgets(str, sizeof(str), fp) || feof(fp)) return 0;
@@ -2918,9 +2920,9 @@ int GModel::readDIFF(const std::string &name)
else
format_read_bi += " %d";
if(sscanf(str, format_read_bi.c_str(), &bi[i]) != 1) return 0;
- Msg::Info("bi[%d]=%d", i, bi[i]);
+ Msg::Info("bi[%d]=%d", i, bi[i]);
}
-
+
while(str[0] != '#'){
if(!fgets(str, sizeof(str), fp) || feof(fp))
break;
@@ -2936,7 +2938,7 @@ int GModel::readDIFF(const std::string &name)
if(!fgets(str, sizeof(str), fp)) return 0;
double xyz[3];
int tmp;
- if(sscanf(str, "%d ( %lf , %lf , %lf ) [%d]", &num, &xyz[0], &xyz[1], &xyz[2],
+ if(sscanf(str, "%d ( %lf , %lf , %lf ) [%d]", &num, &xyz[0], &xyz[1], &xyz[2],
&tmp) != 5) return 0;
elementary[i].resize(tmp + 1);
elementary[i][0] = tmp;
@@ -2946,16 +2948,16 @@ int GModel::readDIFF(const std::string &name)
Msg::Warning("Skipping duplicate vertex %d", num);
else
vertexMap[num] = new MVertex(xyz[0], xyz[1], xyz[2], 0, num);
- if(numVertices > 100000)
+ if(numVertices > 100000)
Msg::ProgressMeter(i + 1, numVertices, "Reading nodes");
// If the vertex numbering is dense, transfer the map into a
// vector to speed up element creation
- if((int)vertexMap.size() == numVertices &&
+ 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)
+ if(minVertex == 1)
vertexVector[0] = 0;
else
vertexVector[numVertices] = 0;
@@ -2964,7 +2966,7 @@ int GModel::readDIFF(const std::string &name)
vertexVector[it->first] = it->second;
vertexMap.clear();
}
- Msg::Info("%d ( %lf , %lf , %lf ) [%d]",i, xyz[0], xyz[1], xyz[2],
+ Msg::Info("%d ( %lf , %lf , %lf ) [%d]",i, xyz[0], xyz[1], xyz[2],
elementary[i][0]);
std::string format_read_bi = "%*d ( %*lf , %*lf , %*lf ) [%*d]";
for(int j = 0; j < elementary[i][0]; j++){
@@ -2974,9 +2976,9 @@ int GModel::readDIFF(const std::string &name)
}
else
format_read_bi += " %d";
- if(sscanf(str, format_read_bi.c_str(), &(elementary[i][j + 1])) != 1)
+ if(sscanf(str, format_read_bi.c_str(), &(elementary[i][j + 1])) != 1)
return 0;
- Msg::Info("elementary[%d][%d]=%d", i + 1, j + 1, elementary[i][j + 1]);
+ Msg::Info("elementary[%d][%d]=%d", i + 1, j + 1, elementary[i][j + 1]);
}
}
while(str[0] != '#'){
@@ -3000,31 +3002,31 @@ int GModel::readDIFF(const std::string &name)
eleType = std::string(eleTypec);
int k2; // local number for the element
int NoVertices; // number of vertices per element
- if(eleType == "ElmT3n2D"){
+ if(eleType == "ElmT3n2D"){
NoVertices = 3;
static int map[3] = {0, 1, 2}; // identical to gmsh
mapping=std::vector<int>(map, map + sizeof(map) / sizeof(int));
type = MSH_TRI_3;
}
- else if(eleType == "ElmT6n2D"){
+ else if(eleType == "ElmT6n2D"){
NoVertices = 6;
static int map[6] = {0, 1, 2, 3, 4, 5}; // identical to gmsh
mapping = std::vector<int>(map, map + sizeof(map) / sizeof(int));
type = MSH_TRI_6;
}
- else if(eleType == "ElmB4n2D"){
+ else if(eleType == "ElmB4n2D"){
NoVertices = 4;
static int map[4] = {0, 1, 3, 2}; // local numbering
mapping = std::vector<int>(map, map + sizeof(map) / sizeof(int));
type = MSH_QUA_4;
}
- else if(eleType == "ElmB8n2D"){
+ else if(eleType == "ElmB8n2D"){
NoVertices = 8;
static int map[8] = {0, 1, 3, 2, 4, 6, 7, 5}; // local numbering
mapping = std::vector<int>(map, map + sizeof(map) / sizeof(int));
type = MSH_QUA_8;
}
- else if(eleType == "ElmB9n2D"){
+ else if(eleType == "ElmB9n2D"){
NoVertices = 9;
static int map[9] = {0, 4, 1, 7, 8, 5, 3, 6, 2}; // local numbering
mapping = std::vector<int>(map, map + sizeof(map) / sizeof(int));
@@ -3035,8 +3037,8 @@ int GModel::readDIFF(const std::string &name)
static int map[4] = {0, 1, 2, 3}; // identical to gmsh
mapping = std::vector<int>(map, map + sizeof(map) / sizeof(int));
type = MSH_TET_4;
- }
- else if(eleType == "ElmT10n3D"){
+ }
+ else if(eleType == "ElmT10n3D"){
NoVertices = 10;
static int map[10] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9}; // local numbering
mapping = std::vector<int>(map, map + sizeof(map) / sizeof(int));
@@ -3047,21 +3049,21 @@ int GModel::readDIFF(const std::string &name)
static int map[8] = {4, 5, 0, 1, 7, 6, 3, 2};
mapping = std::vector<int>(map, map + sizeof(map) / sizeof(int));
type = MSH_HEX_8;
- }
+ }
else if(eleType == "ElmB20n3D"){
NoVertices = 20;
- static int map[20] = {4, 5, 0, 1, 7, 6, 3, 2, 16, 8, 19, 13, 15, 12,
+ static int map[20] = {4, 5, 0, 1, 7, 6, 3, 2, 16, 8, 19, 13, 15, 12,
14, 17, 18, 9, 11};
mapping = std::vector<int>(map, map + sizeof(map) / sizeof(int));
type = MSH_HEX_20;
- }
+ }
else if(eleType == "ElmB27n3D"){
NoVertices = 27;
- static int map[27] = {4, 16, 5, 10, 21, 12, 0, 8, 1, 17, 25, 18, 22,
+ static int map[27] = {4, 16, 5, 10, 21, 12, 0, 8, 1, 17, 25, 18, 22,
26, 23, 9, 20, 11, 7, 19, 6, 15, 24, 14, 3, 13, 2};
mapping = std::vector<int>(map, map + sizeof(map) / sizeof(int));
type = MSH_HEX_27;
- }
+ }
else{
return 0;
}
@@ -3070,11 +3072,11 @@ int GModel::readDIFF(const std::string &name)
if(format_read_vertices[format_read_vertices.size()-2] != '*') {
format_read_vertices[format_read_vertices.size()-1] = '*';
format_read_vertices += "d %d";
- }
+ }
else
format_read_vertices += " %d";
k2 = mapping[k];
- if(sscanf(str, format_read_vertices.c_str(), &ElementsNodes[i-1][k2]) != 1)
+ if(sscanf(str, format_read_vertices.c_str(), &ElementsNodes[i-1][k2]) != 1)
return 0;
}
mapping.clear();
@@ -3082,21 +3084,21 @@ int GModel::readDIFF(const std::string &name)
indices[j] = ElementsNodes[i - 1][j];
std::vector<MVertex*> vertices;
if(vertexVector.size()){
- if(!getVertices(numVerticesPerElement, indices, vertexVector, vertices))
+ if(!getVertices(numVerticesPerElement, indices, vertexVector, vertices))
return 0;
}
else{
- if(!getVertices(numVerticesPerElement, indices, vertexMap, vertices))
+ if(!getVertices(numVerticesPerElement, indices, vertexMap, vertices))
return 0;
}
- createElementMSH(this, num, type, physical, elementary[i-1][1], partition,
- vertices, elements, physicals);
+ createElementMSH(this, num, type, physical, elementary[i-1][1], partition,
+ vertices, elements, physicals);
// trouble if elementary[i-1][0]>1 nodal post-processing needed ?
- if(numElements > 100000)
+ if(numElements > 100000)
Msg::ProgressMeter(i + 1, numElements, "Reading elements");
}
}
-
+
// 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++)
@@ -3132,7 +3134,7 @@ int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll,
Msg::Error("Unable to open file '%s'", name.c_str());
return 0;
}
-
+
if(noPhysicalGroups()) saveAll = true;
// get the number of vertices and index the vertices in a continuous
@@ -3178,7 +3180,7 @@ int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll,
if(entities[i]->physicals.size() || saveAll)
for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++)
dim = std::max(dim, entities[i]->getMeshElement(j)->getDim());
-
+
// loop over all elements we need to save
int numElements = 0, maxNumNodesPerElement = 0;
for(unsigned int i = 0; i < entities.size(); i++){
@@ -3206,7 +3208,7 @@ int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll,
for(std::list<int>::iterator it = boundaryIndicators.begin();
it != boundaryIndicators.end(); it++)
fprintf(fp, " %d", *it);
-
+
fprintf(fp, "\n\n\n");
fprintf(fp," Nodal coordinates and nodal boundary indicators,\n");
fprintf(fp," the columns contain:\n");
@@ -3215,7 +3217,7 @@ int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll,
fprintf(fp," - no of boundary indicators that are set (ON)\n");
fprintf(fp," - the boundary indicators that are set (ON) if any.\n");
fprintf(fp,"#\n");
-
+
// write mesh vertices
for(unsigned int i = 0; i < entities.size(); i++){
for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++){
@@ -3230,7 +3232,7 @@ int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll,
}
}
}
-
+
fprintf(fp, "\n");
fprintf(fp, "\n");
fprintf(fp, " Element types and connectivity\n");
@@ -3253,12 +3255,12 @@ int GModel::writeDIFF(const std::string &name, bool binary, bool saveAll,
}
}
fprintf(fp, "\n");
-
+
fclose(fp);
return 1;
}
-GModel *GModel::createGModel(std::map<int, MVertex*> &vertexMap,
+GModel *GModel::createGModel(std::map<int, MVertex*> &vertexMap,
std::vector<int> &elementNum,
std::vector<std::vector<int> > &vertexIndices,
std::vector<int> &elementType,
@@ -3270,10 +3272,10 @@ GModel *GModel::createGModel(std::map<int, MVertex*> &vertexMap,
std::map<int, std::vector<MElement*> > elements[10];
std::map<int, std::map<int, std::string> > physicals[4];
std::vector<MVertex*> vertexVector;
-
+
int numVertices = (int)vertexMap.size();
int numElement = (int)elementNum.size();
-
+
if(numElement != (int)vertexIndices.size())
Msg::Error("Dimension in vertices numbers");
if(numElement != (int)elementType.size())
@@ -3287,11 +3289,11 @@ GModel *GModel::createGModel(std::map<int, MVertex*> &vertexMap,
std::map<int, MVertex*>::const_iterator it = vertexMap.begin();
std::map<int, MVertex*>::const_iterator end = vertexMap.end();
-
+
int maxVertex = std::numeric_limits<int>::min();
int minVertex = std::numeric_limits<int>::max();
int num;
-
+
for(it; it != end; ++it){
num = it->first;
minVertex = std::min(minVertex, num);
@@ -3299,14 +3301,14 @@ GModel *GModel::createGModel(std::map<int, MVertex*> &vertexMap,
}
if(minVertex == std::numeric_limits<int>::max())
Msg::Error("Could not determine the min index of vertices");
-
+
// If the vertex numbering is dense, transfer the map into a
// vector to speed up element creation
if((minVertex == 1 && maxVertex == numVertices) ||
(minVertex == 0 && maxVertex == numVertices - 1)){
Msg::Info("Vertex numbering is dense");
vertexVector.resize(vertexMap.size() + 1);
- if(minVertex == 1)
+ if(minVertex == 1)
vertexVector[0] = 0;
else
vertexVector[numVertices] = 0;
@@ -3315,7 +3317,7 @@ GModel *GModel::createGModel(std::map<int, MVertex*> &vertexMap,
vertexVector[it->first] = it->second;
vertexMap.clear();
}
-
+
int *indices;
int nbVertices;
for(int i = 0; i < numElement; ++i){
@@ -3337,28 +3339,28 @@ GModel *GModel::createGModel(std::map<int, MVertex*> &vertexMap,
return 0;
}
}
-
+
createElementMSH(gm, num, elementType[i], physical[i], elementary[i],
partition[i], vertices, elements, physicals);
}
-
+
// 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++)
gm->_storeElementsInEntities(elements[i]);
-
+
// associate the correct geometrical entity with each mesh vertex
gm->_associateEntityWithMeshVertices();
-
+
// store the vertices in their associated geometrical entity
if(vertexVector.size())
gm->_storeVerticesInEntities(vertexVector);
else
gm->_storeVerticesInEntities(vertexMap);
-
+
// store the physical tags
for(int i = 0; i < 4; i++)
gm->_storePhysicalTagsInEntities(i, physicals[i]);
-
+
return gm;
}
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 47fba51064cb777acd2f8d20d2e4bc3f0347bc53..082a3841081be467c6fcb43c6677b341dd36045f 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -29,7 +29,7 @@ int MElement::_globalNum = 0;
double MElement::_isInsideTolerance = 1.e-6;
double MElementLessThanLexicographic::tolerance = 1.e-6;
-MElement::MElement(int num, int part) : _visible(1)
+MElement::MElement(int num, int part) : _visible(1)
{
#pragma omp critical
{
@@ -40,12 +40,12 @@ MElement::MElement(int num, int part) : _visible(1)
else{
_num = ++_globalNum;
}
- _partition = (short)part;
+ _partition = (short)part;
}
}
-void MElement::_getEdgeRep(MVertex *v0, MVertex *v1,
- double *x, double *y, double *z, SVector3 *n,
+void MElement::_getEdgeRep(MVertex *v0, MVertex *v1,
+ double *x, double *y, double *z, SVector3 *n,
int faceIndex)
{
x[0] = v0->x(); y[0] = v0->y(); z[0] = v0->z();
@@ -59,10 +59,10 @@ void MElement::_getEdgeRep(MVertex *v0, MVertex *v1,
}
}
-void MElement::_getFaceRep(MVertex *v0, MVertex *v1, MVertex *v2,
+void MElement::_getFaceRep(MVertex *v0, MVertex *v1, MVertex *v2,
double *x, double *y, double *z, SVector3 *n)
{
- x[0] = v0->x(); x[1] = v1->x(); x[2] = v2->x();
+ x[0] = v0->x(); x[1] = v1->x(); x[2] = v2->x();
y[0] = v0->y(); y[1] = v1->y(); y[2] = v2->y();
z[0] = v0->z(); z[1] = v1->z(); z[2] = v2->z();
SVector3 t1(x[1] - x[0], y[1] - y[0], z[1] - z[0]);
@@ -75,7 +75,7 @@ void MElement::_getFaceRep(MVertex *v0, MVertex *v1, MVertex *v2,
char MElement::getVisibility() const
{
if(CTX::instance()->hideUnselected && _visible < 2) return false;
- return _visible;
+ return _visible;
}
double MElement::minEdge()
@@ -115,7 +115,7 @@ void MElement::getShapeFunctions(double u, double v, double w, double s[], int o
else Msg::Error("Function space not implemented for this type of element");
}
-void MElement::getGradShapeFunctions(double u, double v, double w, double s[][3],
+void MElement::getGradShapeFunctions(double u, double v, double w, double s[][3],
int o)
{
const polynomialBasis* fs = getFunctionSpace(o);
@@ -123,6 +123,14 @@ void MElement::getGradShapeFunctions(double u, double v, double w, double s[][3]
else Msg::Error("Function space not implemented for this type of element");
}
+void MElement::getHessShapeFunctions(double u, double v, double w, double s[][3][3],
+ int o)
+{
+ const polynomialBasis* fs = getFunctionSpace(o);
+ if(fs) fs->ddf(u, v, w, s);
+ else Msg::Error("Function space not implemented for this type of element");
+}
+
SPoint3 MElement::barycenter()
{
SPoint3 p(0., 0., 0.);
@@ -149,7 +157,7 @@ std::string MElement::getInfoString()
static double _computeDeterminantAndRegularize(MElement *ele, double jac[3][3])
{
double dJ = 0;
-
+
switch (ele->getDim()) {
case 0:
@@ -157,10 +165,10 @@ static double _computeDeterminantAndRegularize(MElement *ele, double jac[3][3])
dJ = 1.0;
jac[0][0] = jac[1][1] = jac[2][2] = 1.0;
jac[0][1] = jac[1][0] = jac[2][0] = 0.0;
- jac[0][2] = jac[1][2] = jac[2][1] = 0.0;
+ jac[0][2] = jac[1][2] = jac[2][1] = 0.0;
break;
- }
- case 1:
+ }
+ case 1:
{
dJ = sqrt(SQU(jac[0][0]) + SQU(jac[0][1]) + SQU(jac[0][2]));
@@ -179,8 +187,8 @@ static double _computeDeterminantAndRegularize(MElement *ele, double jac[3][3])
norme(b);
prodve(a, b, c);
norme(c);
- jac[0][1] = b[0]; jac[1][1] = b[1]; jac[2][1] = b[2];
- jac[0][2] = c[0]; jac[1][2] = c[1]; jac[2][2] = c[2];
+ jac[0][1] = b[0]; jac[1][1] = b[1]; jac[2][1] = b[2];
+ jac[0][2] = c[0]; jac[1][2] = c[1]; jac[2][2] = c[2];
break;
}
case 2:
@@ -210,7 +218,7 @@ static double _computeDeterminantAndRegularize(MElement *ele, double jac[3][3])
break;
}
}
- return dJ;
+ return dJ;
}
double MElement::getJacobian(double u, double v, double w, double jac[3][3])
@@ -265,7 +273,7 @@ void MElement::pnt(double u, double v, double w, SPoint3 &p)
x += sf[j] * v->x();
y += sf[j] * v->y();
z += sf[j] * v->z();
- }
+ }
p = SPoint3(x, y, z);
}
@@ -306,7 +314,7 @@ void MElement::xyz2uvw(double xyz[3], double uvw[3])
}
double inv[3][3];
inv3x3(jac, inv);
- double un = uvw[0] + inv[0][0] * (xyz[0] - xn) +
+ double un = uvw[0] + inv[0][0] * (xyz[0] - xn) +
inv[1][0] * (xyz[1] - yn) + inv[2][0] * (xyz[2] - zn);
double vn = uvw[1] + inv[0][1] * (xyz[0] - xn) +
inv[1][1] * (xyz[1] - yn) + inv[2][1] * (xyz[2] - zn);
@@ -384,7 +392,7 @@ double MElement::interpolateDiv(double val[], double u, double v, double w,
return fx[0] + fy[1] + fz[2];
}
-void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
+void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
int elementary, int physical, int parentNum,
std::vector<short> *ghosts)
{
@@ -451,8 +459,8 @@ void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
if(physical < 0) revert();
}
-void MElement::writePOS(FILE *fp, bool printElementary, bool printElementNumber,
- bool printGamma, bool printEta, bool printRho,
+void MElement::writePOS(FILE *fp, bool printElementary, bool printElementNumber,
+ bool printGamma, bool printEta, bool printRho,
bool printDisto, double scalingFactor, int elementary)
{
const char *str = getStringForPOS();
@@ -463,7 +471,7 @@ void MElement::writePOS(FILE *fp, bool printElementary, bool printElementNumber,
fprintf(fp, "%s(", str);
for(int i = 0; i < n; i++){
if(i) fprintf(fp, ",");
- fprintf(fp, "%g,%g,%g", getVertex(i)->x() * scalingFactor,
+ fprintf(fp, "%g,%g,%g", getVertex(i)->x() * scalingFactor,
getVertex(i)->y() * scalingFactor, getVertex(i)->z() * scalingFactor);
}
fprintf(fp, "){");
@@ -520,9 +528,9 @@ void MElement::writeSTL(FILE *fp, bool binary, double scalingFactor)
fprintf(fp, "facet normal %g %g %g\n", n[0], n[1], n[2]);
fprintf(fp, " outer loop\n");
for(int j = 0; j < 3; j++)
- fprintf(fp, " vertex %g %g %g\n",
- getVertex(j)->x() * scalingFactor,
- getVertex(j)->y() * scalingFactor,
+ fprintf(fp, " vertex %g %g %g\n",
+ getVertex(j)->x() * scalingFactor,
+ getVertex(j)->y() * scalingFactor,
getVertex(j)->z() * scalingFactor);
fprintf(fp, " endloop\n");
fprintf(fp, "endfacet\n");
@@ -530,9 +538,9 @@ void MElement::writeSTL(FILE *fp, bool binary, double scalingFactor)
fprintf(fp, "facet normal %g %g %g\n", n[0], n[1], n[2]);
fprintf(fp, " outer loop\n");
for(int j = 0; j < 3; j++)
- fprintf(fp, " vertex %g %g %g\n",
- getVertex(qid[j])->x() * scalingFactor,
- getVertex(qid[j])->y() * scalingFactor,
+ fprintf(fp, " vertex %g %g %g\n",
+ getVertex(qid[j])->x() * scalingFactor,
+ getVertex(qid[j])->y() * scalingFactor,
getVertex(qid[j])->z() * scalingFactor);
fprintf(fp, " endloop\n");
fprintf(fp, "endfacet\n");
@@ -622,22 +630,22 @@ void MElement::writeUNV(FILE *fp, int num, int elementary, int physical)
if(physical < 0) revert();
}
-void MElement::writeMESH(FILE *fp, int elementTagType, int elementary,
+void MElement::writeMESH(FILE *fp, int elementTagType, int elementary,
int physical)
{
setVolumePositive();
for(int i = 0; i < getNumVertices(); i++)
fprintf(fp, " %d", getVertex(i)->getIndex());
- fprintf(fp, " %d\n", (elementTagType == 3) ? _partition :
+ fprintf(fp, " %d\n", (elementTagType == 3) ? _partition :
(elementTagType == 2) ? physical : elementary);
}
-void MElement::writeFEA(FILE *fp, int elementTagType, int num, int elementary,
+void MElement::writeFEA(FILE *fp, int elementTagType, int num, int elementary,
int physical)
{
int numVert = getNumVertices();
setVolumePositive();
- fprintf(fp, "%d %d %d", num, (elementTagType == 3) ? _partition :
+ fprintf(fp, "%d %d %d", num, (elementTagType == 3) ? _partition :
(elementTagType == 2) ? physical : elementary, numVert);
for(int i = 0; i < numVert; i++)
fprintf(fp, " %d", getVertex(i)->getIndex());
@@ -654,8 +662,8 @@ void MElement::writeBDF(FILE *fp, int format, int elementTagType, int elementary
int n = getNumVertices();
const char *cont[4] = {"E", "F", "G", "H"};
int ncont = 0;
-
- int tag = (elementTagType == 3) ? _partition : (elementTagType == 2) ?
+
+ int tag = (elementTagType == 3) ? _partition : (elementTagType == 2) ?
physical : elementary;
if(format == 0){ // free field format
@@ -749,14 +757,14 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name)
case MSH_PYR_13 : if(name) *name = "Pyramid 13"; return 5 + 8;
case MSH_PYR_14 : if(name) *name = "Pyramid 14"; return 5 + 8 + 1;
case MSH_POLYH_ : if(name) *name = "Polyhedron"; return 0;
- default:
- Msg::Error("Unknown type of element %d", typeMSH);
- if(name) *name = "Unknown";
+ default:
+ Msg::Error("Unknown type of element %d", typeMSH);
+ if(name) *name = "Unknown";
return 0;
- }
+ }
}
-MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
+MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
int num, int part)
{
switch (type) {
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 5e5af4fc6410ac15caf566111d534b55c7894713..9292dee6071a44253ee83bb445d426ecf163a2fb 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -21,7 +21,7 @@ class GFace;
class binding;
// A mesh element.
-class MElement
+class MElement
{
private:
// the maximum element id number in the mesh
@@ -37,10 +37,10 @@ class MElement
// the tolerance used to determine if a point is inside an element,
// in parametric coordinates
static double _isInsideTolerance;
- void _getEdgeRep(MVertex *v0, MVertex *v1,
+ void _getEdgeRep(MVertex *v0, MVertex *v1,
double *x, double *y, double *z, SVector3 *n,
int faceIndex=-1);
- void _getFaceRep(MVertex *v0, MVertex *v1, MVertex *v2,
+ void _getFaceRep(MVertex *v0, MVertex *v1, MVertex *v2,
double *x, double *y, double *z, SVector3 *n);
public :
MElement(int num=0, int part=0);
@@ -187,13 +187,13 @@ class MElement
virtual double getVolume(){ return 0.; }
virtual int getVolumeSign(){ return 1; }
virtual void setVolumePositive(){ if(getVolumeSign() < 0) revert(); }
-
+
// return an information string for the element
virtual std::string getInfoString();
// get the function space for the element
virtual const polynomialBasis* getFunctionSpace(int order=-1) const { return 0; }
-
+
// return the interpolating nodal shape functions evaluated at point
// (u,v,w) in parametric coordinates (if order == -1, use the
// polynomial order of the element)
@@ -205,17 +205,18 @@ class MElement
// polynomial order of the element)
virtual void getGradShapeFunctions(double u, double v, double w, double s[][3],
int order=-1);
-
+ virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3],
+ int order=-1);
// return the Jacobian of the element evaluated at point (u,v,w) in
// parametric coordinates
double getJacobian(double u, double v, double w, double jac[3][3]);
double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
-
+
// get the point in cartesian coordinates corresponding to the point
// (u,v,w) in parametric coordinates
virtual void pnt(double u, double v, double w, SPoint3 &p);
virtual void primaryPnt(double u, double v, double w, SPoint3 &p);
-
+
// invert the parametrisation
virtual void xyz2uvw(double xyz[3], double uvw[3]);
@@ -225,7 +226,7 @@ class MElement
// interpolate the given nodal data (resp. its gradient, curl and
// divergence) at point (u,v,w) in parametric coordinates
- double interpolate(double val[], double u, double v, double w, int stride=1,
+ double interpolate(double val[], double u, double v, double w, int stride=1,
int order=-1);
void interpolateGrad(double val[], double u, double v, double w, double f[3],
int stride=1, double invjac[3][3]=0, int order=-1);
@@ -244,22 +245,23 @@ class MElement
virtual void writeMSH(FILE *fp, double version=1.0, bool binary=false,
int num=0, int elementary=1, int physical=1,
int parentNum=0, std::vector<short> *ghosts=0);
- virtual void writePOS(FILE *fp, bool printElementary, bool printElementNumber,
- bool printGamma, bool printEta, bool printRho,
+
+ virtual void writePOS(FILE *fp, bool printElementary, bool printElementNumber,
+ bool printGamma, bool printEta, bool printRho,
bool printDisto,double scalingFactor=1.0, int elementary=1);
virtual void writeSTL(FILE *fp, bool binary=false, double scalingFactor=1.0);
virtual void writeVRML(FILE *fp);
virtual void writeUNV(FILE *fp, int num=0, int elementary=1, int physical=1);
virtual void writeVTK(FILE *fp, bool binary=false, bool bigEndian=false);
- virtual void writeMESH(FILE *fp, int elementTagType=1, int elementary=1,
+ virtual void writeMESH(FILE *fp, int elementTagType=1, int elementary=1,
int physical=0);
- virtual void writeFEA(FILE *fp, int elementTagType, int num, int elementary,
+ virtual void writeFEA(FILE *fp, int elementTagType, int num, int elementary,
int physical);
- virtual void writeBDF(FILE *fp, int format=0, int elementTagType=1,
+ virtual void writeBDF(FILE *fp, int format=0, int elementTagType=1,
int elementary=1, int physical=0);
- virtual void writeDIFF(FILE *fp, int num, bool binary=false,
+ virtual void writeDIFF(FILE *fp, int num, bool binary=false,
int physical_property=1);
-
+
// info for specific IO formats (returning 0 means that the element
// is not implemented in that format)
virtual int getTypeForMSH() const { return 0; }
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 08b7c47ca3392a0fbd047f15929818687272b96e..2f53c3098bb32b653e59519d44a645cd8d9e4aec 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -12,13 +12,13 @@
#include "Geo.h"
#include "GFaceCompound.h"
#include "Context.h"
-#include "Os.h"
+#include "OS.h"
discreteFace::discreteFace(GModel *model, int num) : GFace(model, num)
{
Surface *s = Create_Surface(num, MSH_SURF_DISCRETE);
Tree_Add(model->getGEOInternals()->Surfaces, &s);
- meshStatistics.status = GFace::DONE;
+ meshStatistics.status = GFace::DONE;
}
void discreteFace::findEdges(std::map<MEdge, std::vector<int>, Less_Edge> &map_edges)
@@ -36,11 +36,11 @@ void discreteFace::findEdges(std::map<MEdge, std::vector<int>, Less_Edge> &map_e
}
// for the boundary edges, associate the tag of the current discrete face
- for (std::set<MEdge, Less_Edge>::iterator itv = bound_edges.begin();
+ for (std::set<MEdge, Less_Edge>::iterator itv = bound_edges.begin();
itv != bound_edges.end(); ++itv){
std::map<MEdge, std::vector<int>, Less_Edge >::iterator itmap = map_edges.find(*itv);
if (itmap == map_edges.end()){
- std::vector<int> tagFaces;
+ std::vector<int> tagFaces;
tagFaces.push_back(tag());
map_edges.insert(std::make_pair(*itv, tagFaces));
}
@@ -62,7 +62,7 @@ void discreteFace::setBoundEdges(std::vector<int> tagEdges)
}
}
-GPoint discreteFace::point(double par1, double par2) const
+GPoint discreteFace::point(double par1, double par2) const
{
Msg::Error("Cannot evaluate point on discrete face");
return GPoint();
@@ -119,7 +119,7 @@ Pair<SVector3, SVector3> discreteFace::firstDer(const SPoint2 ¶m) const
}
-void discreteFace::secondDer(const SPoint2 ¶m,
+void discreteFace::secondDer(const SPoint2 ¶m,
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
{
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index b2a59a6690968c9da0e3077d7f3d40fc79d2bcba..95dbdb1f1153cc10345a39f4e9560efcc6c19e31 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -53,7 +53,7 @@ static fullMatrix<double> generatePascalSerendipityTriangle(int order)
monomials(0, 0) = 0;
monomials(0, 1) = 0;
-
+
int index = 1;
for (int i = 1; i <= order; i++) {
if (i == order) {
@@ -92,7 +92,7 @@ static fullMatrix<double> generatePascalQuad(int order)
}
}
return monomials;
-}
+}
/*
00 10 20 30 40 ⋯
01 11 21 31 41 ⋯
@@ -157,7 +157,7 @@ static fullMatrix<double> generatePascalQuadSerendip(int order)
static fullMatrix<double> generateMonomialSubspace(int dim, int p)
{
fullMatrix<double> monomials;
-
+
switch (dim) {
case 1:
monomials = fullMatrix<double>(1, 1);
@@ -190,7 +190,7 @@ static fullMatrix<double> generateMonomialSubspace(int dim, int p)
static fullMatrix<double> generatePascalSerendipityTetrahedron(int order)
{
- int nbMonomials = 4 + 6 * std::max(0, order - 1) +
+ int nbMonomials = 4 + 6 * std::max(0, order - 1) +
4 * std::max(0, (order - 2) * (order - 1) / 2);
fullMatrix<double> monomials(nbMonomials, 3);
@@ -231,7 +231,7 @@ static fullMatrix<double> generatePascalTetrahedron(int order)
}
return monomials;
-}
+}
// generate Pascal prism
@@ -254,7 +254,7 @@ static fullMatrix<double> generatePascalPrism(int order)
}
monomials.print("Pri monoms");
return monomials;
-}
+}
static int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
@@ -268,29 +268,29 @@ static void nodepositionface0(int order, double *u, double *v, double *w)
{
int ndofT = nbdoftriangle(order);
if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
-
+
u[0]= 0.; v[0]= 0.; w[0] = 0.;
u[1]= 0.; v[1]= order; w[1] = 0.;
u[2]= order; v[2]= 0.; w[2] = 0.;
// edges
for (int k = 0; k < (order - 1); k++){
- u[3 + k] = 0.;
- v[3 + k] = k + 1;
+ u[3 + k] = 0.;
+ v[3 + k] = k + 1;
w[3 + k] = 0.;
u[3 + order - 1 + k] = k + 1;
- v[3 + order - 1 + k] = order - 1 - k ;
+ v[3 + order - 1 + k] = order - 1 - k ;
w[3 + order - 1 + k] = 0.;
u[3 + 2 * (order - 1) + k] = order - 1 - k;
v[3 + 2 * (order - 1) + k] = 0.;
w[3 + 2 * (order - 1) + k] = 0.;
}
-
+
if (order > 2){
int nbdoftemp = nbdoftriangle(order - 3);
- nodepositionface0(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
+ nodepositionface0(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
&w[3 + 3* (order - 1)]);
for (int k = 0; k < nbdoftemp; k++){
u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
@@ -310,23 +310,23 @@ static void nodepositionface1(int order, double *u, double *v, double *w)
{
int ndofT = nbdoftriangle(order);
if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
-
+
u[0] = 0.; v[0]= 0.; w[0] = 0.;
u[1] = order; v[1]= 0.; w[1] = 0.;
u[2] = 0.; v[2]= 0.; w[2] = order;
// edges
for (int k = 0; k < (order - 1); k++){
- u[3 + k] = k + 1;
- v[3 + k] = 0.;
+ u[3 + k] = k + 1;
+ v[3 + k] = 0.;
w[3 + k] = 0.;
-
+
u[3 + order - 1 + k] = order - 1 - k;
v[3 + order - 1 + k] = 0.;
- w[3 + order - 1+ k ] = k + 1;
-
+ w[3 + order - 1+ k ] = k + 1;
+
u[3 + 2 * (order - 1) + k] = 0. ;
v[3 + 2 * (order - 1) + k] = 0.;
- w[3 + 2 * (order - 1) + k] = order - 1 - k;
+ w[3 + 2 * (order - 1) + k] = order - 1 - k;
}
if (order > 2){
int nbdoftemp = nbdoftriangle(order - 3);
@@ -350,13 +350,13 @@ static void nodepositionface2(int order, double *u, double *v, double *w)
{
int ndofT = nbdoftriangle(order);
if (order == 0) { u[0] = 0.; v[0] = 0.; return; }
-
+
u[0]= 0.; v[0]= 0.; w[0] = 0.;
u[1]= 0.; v[1]= 0.; w[1] = order;
u[2]= 0.; v[2]= order; w[2] = 0.;
// edges
for (int k = 0; k < (order - 1); k++){
-
+
u[3 + k] = 0.;
v[3 + k] = 0.;
w[3 + k] = k + 1;
@@ -364,7 +364,7 @@ static void nodepositionface2(int order, double *u, double *v, double *w)
u[3 + order - 1 + k] = 0.;
v[3 + order - 1 + k] = k + 1;
w[3 + order - 1 + k] = order - 1 - k;
-
+
u[3 + 2 * (order - 1) + k] = 0.;
v[3 + 2 * (order - 1) + k] = order - 1 - k;
w[3 + 2 * (order - 1) + k] = 0.;
@@ -391,7 +391,7 @@ static void nodepositionface3(int order, double *u, double *v, double *w)
{
int ndofT = nbdoftriangle(order);
if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
-
+
u[0]= 0.; v[0]= 0.; w[0] = order;
u[1]= order; v[1]= 0.; w[1] = 0.;
u[2]= 0.; v[2]= order; w[2] = 0.;
@@ -405,9 +405,9 @@ static void nodepositionface3(int order, double *u, double *v, double *w)
u[3 + order - 1 + k] = order - 1 - k;
v[3 + order - 1 + k] = k + 1;
w[3 + order - 1 + k] = 0.;
-
+
u[3 + 2 * (order - 1) + k] = 0.;
- v[3 + 2 * (order - 1) + k] = order - 1 - k;
+ v[3 + 2 * (order - 1) + k] = order - 1 - k;
w[3 + 2 * (order - 1) + k] = k + 1;
}
if (order > 2){
@@ -427,13 +427,13 @@ static void nodepositionface3(int order, double *u, double *v, double *w)
}
}
-static fullMatrix<double> gmshGeneratePointsTetrahedron(int order, bool serendip)
+static fullMatrix<double> gmshGeneratePointsTetrahedron(int order, bool serendip)
{
- int nbPoints =
+ int nbPoints =
(serendip ?
4 + 6 * std::max(0, order - 1) + 4 * std::max(0, (order - 2) * (order - 1) / 2) :
(order + 1) * (order + 2) * (order + 3) / 6);
-
+
fullMatrix<double> point(nbPoints, 3);
double overOrder = (order == 0 ? 1. : 1. / order);
@@ -446,90 +446,90 @@ static fullMatrix<double> gmshGeneratePointsTetrahedron(int order, bool serendip
point(1, 0) = order;
point(1, 1) = 0;
point(1, 2) = 0;
-
+
point(2, 0) = 0.;
point(2, 1) = order;
point(2, 2) = 0.;
-
+
point(3, 0) = 0.;
point(3, 1) = 0.;
point(3, 2) = order;
// edges e5 and e6 switched in original version, opposite direction
// the template has been defined in table edges_tetra and faces_tetra (MElement.h)
-
+
if (order > 1) {
for (int k = 0; k < (order - 1); k++) {
point(4 + k, 0) = k + 1;
point(4 + order - 1 + k, 0) = order - 1 - k;
- point(4 + 2 * (order - 1) + k, 0) = 0.;
+ point(4 + 2 * (order - 1) + k, 0) = 0.;
point(4 + 3 * (order - 1) + k, 0) = 0.;
// point(4 + 4 * (order - 1) + k, 0) = order - 1 - k;
// point(4 + 5 * (order - 1) + k, 0) = 0.;
point(4 + 4 * (order - 1) + k, 0) = 0.;
point(4 + 5 * (order - 1) + k, 0) = k+1;
-
+
point(4 + k, 1) = 0.;
point(4 + order - 1 + k, 1) = k + 1;
- point(4 + 2 * (order - 1) + k, 1) = order - 1 - k;
+ point(4 + 2 * (order - 1) + k, 1) = order - 1 - k;
point(4 + 3 * (order - 1) + k, 1) = 0.;
// point(4 + 4 * (order - 1) + k, 1) = 0.;
// point(4 + 5 * (order - 1) + k, 1) = order - 1 - k;
point(4 + 4 * (order - 1) + k, 1) = k+1;
point(4 + 5 * (order - 1) + k, 1) = 0.;
-
+
point(4 + k, 2) = 0.;
point(4 + order - 1 + k, 2) = 0.;
- point(4 + 2 * (order - 1) + k, 2) = 0.;
+ point(4 + 2 * (order - 1) + k, 2) = 0.;
point(4 + 3 * (order - 1) + k, 2) = order - 1 - k;
point(4 + 4 * (order - 1) + k, 2) = order - 1 - k;
point(4 + 5 * (order - 1) + k, 2) = order - 1 - k;
}
-
+
if (order > 2) {
int ns = 4 + 6 * (order - 1);
int nbdofface = nbdoftriangle(order - 3);
-
+
double *u = new double[nbdofface];
double *v = new double[nbdofface];
double *w = new double[nbdofface];
-
+
nodepositionface0(order - 3, u, v, w);
// u-v plane
-
+
for (int i = 0; i < nbdofface; i++){
point(ns + i, 0) = u[i] * (order - 3) + 1.;
point(ns + i, 1) = v[i] * (order - 3) + 1.;
point(ns + i, 2) = w[i] * (order - 3);
}
-
+
ns = ns + nbdofface;
// u-w plane
-
+
nodepositionface1(order - 3, u, v, w);
-
+
for (int i=0; i < nbdofface; i++){
point(ns + i, 0) = u[i] * (order - 3) + 1.;
point(ns + i, 1) = v[i] * (order - 3) ;
point(ns + i, 2) = w[i] * (order - 3) + 1.;
}
- // v-w plane
-
+ // v-w plane
+
ns = ns + nbdofface;
nodepositionface2(order - 3, u, v, w);
-
+
for (int i = 0; i < nbdofface; i++){
point(ns + i, 0) = u[i] * (order - 3);
point(ns + i, 1) = v[i] * (order - 3) + 1.;
point(ns + i, 2) = w[i] * (order - 3) + 1.;
}
- // u-v-w plane
-
+ // u-v-w plane
+
ns = ns + nbdofface;
nodepositionface3(order - 3, u, v, w);
@@ -545,9 +545,9 @@ static fullMatrix<double> gmshGeneratePointsTetrahedron(int order, bool serendip
delete [] u;
delete [] v;
delete [] w;
-
+
if (!serendip && order > 3) {
-
+
fullMatrix<double> interior = gmshGeneratePointsTetrahedron(order - 4, false);
for (int k = 0; k < interior.size1(); k++) {
point(ns + k, 0) = 1. + interior(k, 0) * (order - 4);
@@ -558,19 +558,19 @@ static fullMatrix<double> gmshGeneratePointsTetrahedron(int order, bool serendip
}
}
}
-
- point.scale(overOrder);
+
+ point.scale(overOrder);
return point;
}
-static fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip)
+static fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip)
{
int nbPoints = serendip ? 3 * order : (order + 1) * (order + 2) / 2;
fullMatrix<double> point(nbPoints, 2);
-
+
point(0, 0) = 0;
point(0, 1) = 0;
-
+
double dd = 1. / order;
if (order > 0) {
@@ -578,29 +578,29 @@ static fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip)
point(1, 1) = 0;
point(2, 0) = 0;
point(2, 1) = 1;
-
+
int index = 3;
-
+
if (order > 1) {
-
+
double ksi = 0;
double eta = 0;
-
+
for (int i = 0; i < order - 1; i++, index++) {
ksi += dd;
point(index, 0) = ksi;
point(index, 1) = eta;
}
-
+
ksi = 1.;
-
+
for (int i = 0; i < order - 1; i++, index++) {
ksi -= dd;
eta += dd;
point(index, 0) = ksi;
point(index, 1) = eta;
- }
-
+ }
+
eta = 1.;
ksi = 0.;
@@ -608,7 +608,7 @@ static fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip)
eta -= dd;
point(index, 0) = ksi;
point(index, 1) = eta;
- }
+ }
if (order > 2 && !serendip) {
fullMatrix<double> inner = gmshGeneratePointsTriangle(order - 3, serendip);
@@ -618,10 +618,10 @@ static fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip)
}
}
}
- return point;
+ return point;
}
-static fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip)
+static fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip)
{
const double prism18Pts[18][3] = {
{0, 0, -1}, // 0
@@ -644,9 +644,9 @@ static fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip)
{0.5, 0.5, 0}, // 17
};
- int nbPoints = (order + 1)*(order + 1)*(order + 2)/2;
+ int nbPoints = (order + 1)*(order + 1)*(order + 2)/2;
fullMatrix<double> point(nbPoints, 3);
-
+
int index = 0;
fullMatrix<double> triPoints = gmshGeneratePointsTriangle(order,false);
fullMatrix<double> linePoints = generate1DPoints(order);
@@ -655,7 +655,7 @@ static fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip)
for (int i =0; i<18; i++)
for (int j=0; j<3;j++)
point(i,j) = prism18Pts[i][j];
- else
+ else
for (int j = 0; j <linePoints.size1() ; j++) {
for (int i = 0; i < triPoints.size1(); i++) {
point(index,0) = triPoints(i,0);
@@ -664,17 +664,17 @@ static fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip)
index ++;
}
}
-
+
// point.print("Pri ipts");
return point;
}
-static fullMatrix<double> gmshGeneratePointsQuad(int order, bool serendip)
+static fullMatrix<double> gmshGeneratePointsQuad(int order, bool serendip)
{
int nbPoints = serendip ? order*4 : (order+1)*(order+1);
fullMatrix<double> point(nbPoints, 2);
-
+
double dd = 1. / order;
if (order > 0) {
@@ -686,7 +686,7 @@ static fullMatrix<double> gmshGeneratePointsQuad(int order, bool serendip)
point(2, 1) = 1;
point(3, 0) = -1;
point(3, 1) = 1;
-
+
if (order > 1) {
int index = 4;
const static int edges[4][2]={{0,1},{1,2},{2,3},{3,0}};
@@ -704,25 +704,25 @@ static fullMatrix<double> gmshGeneratePointsQuad(int order, bool serendip)
point.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
}
}
- }
+ }
else {
point(0, 0) = 0;
point(0, 1) = 0;
}
- return point;
+ return point;
}
static fullMatrix<double> generateLagrangeMonomialCoefficients
- (const fullMatrix<double>& monomial, const fullMatrix<double>& point)
+ (const fullMatrix<double>& monomial, const fullMatrix<double>& point)
{
if(monomial.size1() != point.size1() || monomial.size2() != point.size2()){
Msg::Fatal("Wrong sizes for Lagrange coefficients generation");
return fullMatrix<double>(1, 1);
}
-
+
int ndofs = monomial.size1();
int dim = monomial.size2();
-
+
fullMatrix<double> Vandermonde(ndofs, ndofs);
for (int i = 0; i < ndofs; i++) {
for (int j = 0; j < ndofs; j++) {
@@ -755,7 +755,7 @@ static void getFaceClosure(int iFace, int iSign, int iRotate, std::vector<int> &
closure[i] = order1node[iFace][k];
}
for (int i = 0;i < 3; ++i){
- int edgenumber = iSign *
+ int edgenumber = iSign *
face[iFace][(6 + i * iSign + (-1 + iSign) / 2 + iRotate) % 3]; //- iSign * iRotate
for (int k = 0; k < (order - 1); k++){
if (edgenumber > 0)
@@ -763,7 +763,7 @@ static void getFaceClosure(int iFace, int iSign, int iRotate, std::vector<int> &
4 + (edgenumber - 1) * (order - 1) + k;
else
closure[3 + i * (order - 1) + k] =
- 4 + (-edgenumber) * (order - 1) - 1 - k;
+ 4 + (-edgenumber) * (order - 1) - 1 - k;
}
}
int fi = 3 + 3 * (order - 1);
@@ -781,26 +781,26 @@ static void getFaceClosure(int iFace, int iSign, int iRotate, std::vector<int> &
for (int l = 0; l < orderint - 1; l++){
for (int ei = 0; ei < 3; ei++){
int edgenumber = (6 + ei * iSign + (-1 + iSign) / 2 + iRotate) % 3; //- iSign * iRotate
- if (iSign > 0)
- closure[fi + ei * (orderint - 1) + l] =
+ if (iSign > 0)
+ closure[fi + ei * (orderint - 1) + l] =
ti + edgenumber * (orderint - 1) + l;
else
closure[fi + ei * (orderint - 1) + l] =
- ti + (1 + edgenumber) * (orderint - 1) - 1 - l;
+ ti + (1 + edgenumber) * (orderint - 1) - 1 - l;
}
}
- fi = fi + 3 * (orderint - 1); ti = ti + 3 * (orderint - 1);
+ fi = fi + 3 * (orderint - 1); ti = ti + 3 * (orderint - 1);
}
else {
closure[fi] = ti;
- ti++;
+ ti++;
fi++;
- }
+ }
}
break;
}
-}
+}
static void generate3dFaceClosure(polynomialBasis::clCont &closure, int order)
{
@@ -810,7 +810,7 @@ static void generate3dFaceClosure(polynomialBasis::clCont &closure, int order)
for (int iSign = 1; iSign >= -1; iSign -= 2){
for (int iFace = 0; iFace < 4; iFace++){
std::vector<int> closure_face;
- getFaceClosure(iFace, iSign, iRotate, closure_face, order);
+ getFaceClosure(iFace, iSign, iRotate, closure_face, order);
closure.push_back(closure_face);
}
}
@@ -856,7 +856,7 @@ static void generate3dFaceClosurePrism(polynomialBasis::clCont &closure, int ord
for (int iSign = 1; iSign >= -1; iSign -= 2){
for (int iFace = 0; iFace < 5; iFace++){
std::vector<int> closure_face;
- getFaceClosurePrism(iFace, iSign, iRotate, closure_face, order);
+ getFaceClosurePrism(iFace, iSign, iRotate, closure_face, order);
closure.push_back(closure_face);
}
}
@@ -887,17 +887,17 @@ static void generate1dVertexClosure(polynomialBasis::clCont &closure)
closure.resize(2);
closure[0].push_back(0);
closure[1].push_back(1);
-
+
}
std::map<int, polynomialBasis> polynomialBases::fs;
-const polynomialBasis &polynomialBases::find(int tag)
+const polynomialBasis &polynomialBases::find(int tag)
{
std::map<int, polynomialBasis>::const_iterator it = fs.find(tag);
if (it != fs.end()) return it->second;
polynomialBasis F;
F.numFaces = -1;
-
+
switch (tag){
case MSH_PNT:
F.numFaces = 1;
@@ -933,7 +933,7 @@ const polynomialBasis &polynomialBases::find(int tag)
F.monomials = generate1DMonomials(5);
F.points = generate1DPoints(5);
generate1dVertexClosure(F.closures);
- break;
+ break;
case MSH_TRI_3 :
F.numFaces = 3;
F.monomials = generatePascalTriangle(1);
@@ -1090,7 +1090,7 @@ const polynomialBasis &polynomialBases::find(int tag)
F.points = gmshGeneratePointsPrism(2, false);
generate3dFaceClosurePrism(F.closures, 2);
break;
-
+
default :
Msg::Error("Unknown function space %d: reverting to TET_4", tag);
F.numFaces = 4;
@@ -1098,7 +1098,7 @@ const polynomialBasis &polynomialBases::find(int tag)
F.points = gmshGeneratePointsTetrahedron(1, false);
generate3dFaceClosure(F.closures, 1);
break;
- }
+ }
F.coefficients = generateLagrangeMonomialCoefficients(F.monomials, F.points);
// printf("Case: %d coeffs:\n",tag);
// for (int i = 0; i<F.coefficients.size1(); i++) {
@@ -1107,7 +1107,7 @@ const polynomialBasis &polynomialBases::find(int tag)
// }
// printf("\n");
// }
-
+
fs.insert(std::make_pair(tag, F));
return fs[tag];
}
@@ -1125,9 +1125,9 @@ const fullMatrix<double> &polynomialBases::findInjector(int tag1, int tag2)
const polynomialBasis& fs2 = find(tag2);
fullMatrix<double> inj(fs1.points.size1(), fs2.points.size1());
-
+
double sf[256];
-
+
for (int i = 0; i < fs1.points.size1(); i++) {
fs2.f(fs1.points(i, 0), fs1.points(i, 1), fs1.points(i, 2), sf);
for (int j = 0; j < fs2.points.size1(); j++) inj(i, j) = sf[j];
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 6055d1b32b0c3774d91ed8d5b9b466397a2f5f28..8ab2429d8ebb8b145d4080b406e77bf1428aecd9 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -10,6 +10,7 @@
#include <map>
#include <vector>
#include "fullMatrix.h"
+#include <iostream>
// presently those function spaces are only for simplices and quads;
// should be extended to other elements like hexes
@@ -31,7 +32,7 @@ class polynomialBasis
inline int getClosureId(int iEl, int iSign=1, int iRot=0) const {
return iEl + numFaces*(iSign == 1 ? 0 : 1) + 2*numFaces*iRot;
}
- inline void evaluateMonomials(double u, double v, double w, double p[]) const
+ inline void evaluateMonomials(double u, double v, double w, double p[]) const
{
for (int j = 0; j < monomials.size1(); j++) {
p[j] = pow(u, (int)monomials(j, 0));
@@ -60,7 +61,7 @@ class polynomialBasis
grads[i][2] = 0;
for(int j = 0; j < coefficients.size2(); j++){
if ((monomials)(j, 0) > 0)
- grads[i][0] += (coefficients)(i, j) *
+ grads[i][0] += (coefficients)(i, j) *
pow(u, (monomials)(j, 0) - 1) * (monomials)(j, 0);
}
}
@@ -108,6 +109,88 @@ class polynomialBasis
break;
}
}
+ inline void ddf(double u, double v, double w, double hess[][3][3]) const
+ {
+ switch (monomials.size2()) {
+ case 1:
+ for (int i = 0; i < coefficients.size1(); i++){
+ hess[i][0][0] = hess[i][0][1] = hess[i][0][2] = 0;
+ hess[i][1][0] = hess[i][1][1] = hess[i][1][2] = 0;
+ hess[i][2][0] = hess[i][2][1] = hess[i][2][2] = 0;
+
+ for(int j = 0; j < coefficients.size2(); j++){
+ if ((monomials)(j, 0) > 1) // second derivative !=0
+ hess[i][0][0] += (coefficients)(i, j) * pow(u, (monomials)(j, 0) - 2) * (monomials)(j, 0) * ((monomials)(j, 0)-1);
+ }
+ }
+ break;
+ case 2:
+ for (int i = 0; i < coefficients.size1(); i++){
+ hess[i][0][0] = hess[i][0][1] = hess[i][0][2] = 0;
+ hess[i][1][0] = hess[i][1][1] = hess[i][1][2] = 0;
+ hess[i][2][0] = hess[i][2][1] = hess[i][2][2] = 0;
+ for(int j = 0; j < coefficients.size2(); j++){
+ if ((monomials)(j, 0) > 1) // second derivative !=0
+ hess[i][0][0] += (coefficients)(i, j) *
+ pow(u, (monomials)(j, 0) - 2) * (monomials)(j, 0) * ((monomials)(j, 0)-1) *
+ pow(v, (monomials)(j, 1));
+ if (((monomials)(j,1)>0) and ((monomials)(j,0)>0))
+ hess[i][0][1]+=(coefficients)(i, j) *
+ pow(u, (monomials)(j, 0) - 1) * (monomials)(j, 0) *
+ pow(v, (monomials)(j, 1)-1) * (monomials)(j, 1);
+ if ((monomials)(j, 1) > 1)
+ hess[i][1][1] += (coefficients)(i, j) *
+ pow(u, (monomials)(j, 0)) *
+ pow(v, (monomials)(j, 1) - 2) * (monomials)(j, 1) * ((monomials)(j, 1)-1);
+ }
+ hess[i][1][0]=hess[i][0][1];
+ }
+ break;
+ case 3:
+ for (int i = 0; i < coefficients.size1(); i++){
+ hess[i][0][0] = hess[i][0][1] = hess[i][0][2] = 0;
+ hess[i][1][0] = hess[i][1][1] = hess[i][1][2] = 0;
+ hess[i][2][0] = hess[i][2][1] = hess[i][2][2] = 0;
+ for(int j = 0; j < coefficients.size2(); j++){
+ if ((monomials)(j, 0) > 1)
+ hess[i][0][0] += (coefficients)(i, j) *
+ pow(u, (monomials)(j, 0) - 2) * (monomials)(j, 0) * ((monomials)(j, 0)-1) *
+ pow(v, (monomials)(j, 1)) *
+ pow(w, (monomials)(j, 2));
+
+ if (((monomials)(j,0)>0) and ((monomials)(j,1)>0))
+ hess[i][0][1] += (coefficients)(i, j) *
+ pow(u, (monomials)(j, 0) - 1) * (monomials)(j, 0) *
+ pow(v, (monomials)(j, 1) - 1) * (monomials)(j, 1) *
+ pow(w, (monomials)(j, 2));
+ if (((monomials)(j,0)>0) and ((monomials)(j,2)>0))
+ hess[i][0][2] += (coefficients)(i, j) *
+ pow(u, (monomials)(j, 0) - 1) * (monomials)(j, 0) *
+ pow(v, (monomials)(j, 1)) *
+ pow(w, (monomials)(j, 2) - 1) * (monomials)(j, 2);
+ if ((monomials)(j, 1) > 1)
+ hess[i][1][1] += (coefficients)(i, j) *
+ pow(u, (monomials)(j, 0)) *
+ pow(v, (monomials)(j, 1) - 2) * (monomials)(j, 1) * ((monomials)(j, 1)-1) *
+ pow(w, (monomials)(j, 2));
+ if (((monomials)(j,1)>0) and ((monomials)(j,2)>0))
+ hess[i][1][2] += (coefficients)(i, j) *
+ pow(u, (monomials)(j, 0)) *
+ pow(v, (monomials)(j, 1) - 1) * (monomials)(j, 1) *
+ pow(w, (monomials)(j, 2) - 1) * (monomials)(j, 2);
+ if ((monomials)(j, 2) > 1)
+ hess[i][2][2] += (coefficients)(i, j) *
+ pow(u, (monomials)(j, 0)) *
+ pow(v, (monomials)(j, 1)) *
+ pow(w, (monomials)(j, 2) - 2) * (monomials)(j, 2) * ((monomials)(j, 2)-1);
+ }
+ hess[i][1][0]=hess[i][0][1];
+ hess[i][2][0]=hess[i][0][2];
+ hess[i][2][1]=hess[i][1][2];
+ }
+ break;
+ }
+ }
};
class polynomialBases
diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index 41f39ef10479d75f43a92f6d64d3657d9c601ce6..6d878a2434caa7c1f9543a142270ae4b9b00a6c2 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -68,7 +68,7 @@ class DofAffineConstraint{
std::vector<std::pair<Dof, typename dofTraits<T>::MatType> > linear;
typename dofTraits<T>::VecType shift;
};
-
+
// A manager for degrees of freedoms, templated on the value of a dof
// (what the functional returns): float, double, complex<double>,
// fullVecor<double>, ...
@@ -82,24 +82,24 @@ class dofManager{
// equations:
// dataMat * DofVec = \sum_i dataMat_i * DofVec_i + dataVec
std::map<std::pair<Dof, dataMat>, DofAffineConstraint<dataVec> > constraints;
-
+
// fixations on full blocks, treated by eliminating equations:
- // DofVec = dataVec
+ // DofVec = dataVec
std::map<Dof, dataVec> fixed;
// initial conditions
std::map<Dof, std::vector<dataVec> > initial;
-
+
// numbering of unknown dof blocks
std::map<Dof, int> unknown;
-
+
// associatations
std::map<Dof, Dof> associatedWith;
-
+
// linearSystems
std::map<const std::string, linearSystem<dataMat>*> _linearSystems;
linearSystem<dataMat> *_current;
-
+
public:
dofManager(linearSystem<dataMat> *l) : _current(l) { _linearSystems["A"] = l; }
dofManager(linearSystem<dataMat> *l1, linearSystem<dataMat> *l2) : _current(l1) {
@@ -290,7 +290,7 @@ class dofManager{
assemble(vR->getNum(), Dof::createTypeWithTwoInts(iCompR, iFieldR),
vC->getNum(), Dof::createTypeWithTwoInts(iCompC, iFieldC),
value);
- }
+ }
inline void assemble(const Dof &R, const dataMat &value)
{
if(!_current->isAllocated()) _current->allocate(unknown.size());
@@ -307,15 +307,15 @@ class dofManager{
const dataMat &value)
{
assemble(vR->getNum(), Dof::createTypeWithTwoInts(iCompR, iFieldR), value);
- }
+ }
int sizeOfR() const { return unknown.size(); }
int sizeOfF() const { return fixed.size(); }
- void systemSolve(){ _current->systemSolve(); }
+ void systemSolve(){ _current->systemSolve(); }
void systemClear()
{
_current->zeroMatrix();
_current->zeroRightHandSide();
- }
+ }
inline void setCurrentMatrix(std::string name){
typename std::map<const std::string, linearSystem<dataMat>*>::iterator it = _linearSystems.find(name);
if(it != _linearSystems.end())
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index 693e6817e609e7ea4fbcf202588086ed9a58a159..a9ba9b975d90c18a77d463988b89fe7629d9248f 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -99,7 +99,7 @@ void elasticitySolver::readInputFile(const std::string &fn)
neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
neu._tag=node;
neu.onWhat=BoundaryCondition::ON_VERTEX;
- allNeumann.push_back(neu);
+ allNeumann.push_back(neu);
}
else if (!strcmp(what, "EdgeForce")){
double val1, val2, val3;
@@ -110,7 +110,7 @@ void elasticitySolver::readInputFile(const std::string &fn)
neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
neu._tag=edge;
neu.onWhat=BoundaryCondition::ON_EDGE;
- allNeumann.push_back(neu);
+ allNeumann.push_back(neu);
}
else if (!strcmp(what, "FaceForce")){
double val1, val2, val3;
@@ -121,7 +121,7 @@ void elasticitySolver::readInputFile(const std::string &fn)
neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
neu._tag=face;
neu.onWhat=BoundaryCondition::ON_FACE;
- allNeumann.push_back(neu);
+ allNeumann.push_back(neu);
}
else if (!strcmp(what, "VolumeForce")){
double val1, val2, val3;
@@ -132,7 +132,7 @@ void elasticitySolver::readInputFile(const std::string &fn)
neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
neu._tag=volume;
neu.onWhat=BoundaryCondition::ON_VOLUME;
- allNeumann.push_back(neu);
+ allNeumann.push_back(neu);
}
else if (!strcmp(what, "MeshFile")){
char name[245];
@@ -160,7 +160,7 @@ void elasticitySolver::solve()
if (pAssembler) delete pAssembler;
pAssembler = new dofManager<double>(lsys);
- // we first do all fixations. the behavior of the dofManager is to
+ // 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
@@ -217,8 +217,8 @@ void elasticitySolver::solve()
#if defined(HAVE_POST)
-static double vonMises(dofManager<double> *a, MElement *e,
- double u, double v, double w,
+static double vonMises(dofManager<double> *a, MElement *e,
+ double u, double v, double w,
double E, double nu, int _tag)
{
double valx[256];
@@ -279,7 +279,7 @@ PView* elasticitySolver::buildDisplacementView (const std::string &postFileName)
data[(*it)->getNum()]=vec;
}
PView *pv = new PView (postFileName, "NodeData", pModel, data, 0.0);
- return pv;
+ return pv;
}
PView *elasticitySolver::buildElasticEnergyView(const std::string &postFileName)
diff --git a/Solver/functionSpace.h b/Solver/functionSpace.h
index b25d0e31f294f6e92d1143f3081be0c2dec6eb2f..91fa27586d665bcb57ce30e19e2b9db84d70084c 100644
--- a/Solver/functionSpace.h
+++ b/Solver/functionSpace.h
@@ -25,8 +25,8 @@ template<> struct TensorialTraits<double>
{
typedef double ValType;
typedef SVector3 GradType;
-/* typedef STensor3 HessType;
- typedef SVoid DivType;
+ typedef STensor3 HessType;
+/* typedef SVoid DivType;
typedef SVoid CurlType;*/
};
@@ -34,9 +34,9 @@ template<> struct TensorialTraits<SVector3>
{
typedef SVector3 ValType;
typedef STensor3 GradType;
-/* typedef SVoid HessType;
+ typedef STensor3 HessType;
typedef double DivType;
- typedef SVector3 CurlType;*/
+ typedef SVector3 CurlType;
};
class FunctionSpaceBase
@@ -52,9 +52,12 @@ class FunctionSpace : public FunctionSpaceBase
public:
typedef typename TensorialTraits<T>::ValType ValType;
typedef typename TensorialTraits<T>::GradType GradType;
+ typedef typename TensorialTraits<T>::HessType HessType;
virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)=0;
virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads) =0;
virtual void gradfuvw(MElement *ele, double u, double v, double w,std::vector<GradType> &grads) {}; // should return to pure virtual once all is done.
+ virtual void hessfuvw(MElement *ele, double u, double v, double
+w,std::vector<HessType> &hess)=0;
virtual int getNumKeys(MElement *ele)=0; // if one needs the number of dofs
virtual void getKeys(MElement *ele, std::vector<Dof> &keys)=0;
};
@@ -64,6 +67,7 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
public:
typedef TensorialTraits<double>::ValType ValType;
typedef TensorialTraits<double>::GradType GradType;
+ typedef TensorialTraits<double>::HessType HessType;
protected:
int _iField; // field number (used to build dof keys)
@@ -86,6 +90,8 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
ele->getShapeFunctions(u, v, w, &(vals[curpos]));
}
+
+ // Fonction renvoyant un vecteur contenant le grandient de chaque FF
virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
{
if (ele->getParent()) ele = ele->getParent();
@@ -103,7 +109,22 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2],
invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2]
));
- }
+ };
+ // Fonction renvoyant un vecteur contenant le hessien [][] de chaque FF dans l'espace ISOPARAMETRIQUE
+ virtual void hessfuvw(MElement *ele, double u, double v, double w,std::vector<HessType> &hess)
+ {
+ int ndofs= ele->getNumVertices(); // ATTENTION RETOURNE LE NBBRE DE NOEUDS ET PAS LE NBRE DE DDL
+ hess.reserve(hess.size()+ndofs); // permet de mettre les composantes suivantes à la suite du vecteur
+ double hessuvw[256][3][3];
+ ele->getHessShapeFunctions(u, v, w, hessuvw);
+ HessType hesst;
+ for (int i=0;i<ndofs;++i){
+ hesst(0,0) = hessuvw[i][0][0] ; hesst(0,1) = hessuvw[i][0][1] ; hesst(0,2) = hessuvw[i][0][2];
+ hesst(1,0) = hessuvw[i][1][0] ; hesst(1,1) = hessuvw[i][1][1] ; hesst(1,2) = hessuvw[i][1][2];
+ hesst(2,0) = hessuvw[i][2][0] ; hesst(2,1) = hessuvw[i][2][1] ; hesst(2,2) = hessuvw[i][2][2];
+ hess.push_back(hesst);
+ }
+ };
virtual void gradfuvw(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
{
@@ -138,6 +159,7 @@ template <class T> class ScalarToAnyFunctionSpace : public FunctionSpace<T>
public :
typedef typename TensorialTraits<T>::ValType ValType;
typedef typename TensorialTraits<T>::GradType GradType;
+ typedef typename TensorialTraits<T>::HessType HessType;
protected :
std::vector<T> multipliers;
std::vector<int> comp;
@@ -195,8 +217,10 @@ public :
}
}
}
-
-
+ virtual void hessfuvw(MElement *ele, double u, double v, double w,std::vector<HessType> &hess)
+ {
+ ScalarFS->hessfuvw(ele,u,v,w,hess);
+ }
virtual void gradfuvw(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
{
std::vector<SVector3> gradsd;
@@ -274,6 +298,7 @@ class CompositeFunctionSpace : public FunctionSpace<T>
public:
typedef typename TensorialTraits<T>::ValType ValType;
typedef typename TensorialTraits<T>::GradType GradType;
+ typedef typename TensorialTraits<T>::HessType HessType;
typedef typename std::vector<FunctionSpace<T>* >::iterator iterFS;
protected:
@@ -313,6 +338,12 @@ class CompositeFunctionSpace : public FunctionSpace<T>
(*it)->gradf(ele,u,v,w,grads);
}
+ virtual void hessfuvw(MElement *ele, double u, double v, double w,std::vector<HessType> &hess)
+ {
+ for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
+ (*it)->hessfuvw(ele,u,v,w,hess);
+ }
+
virtual int getNumKeys(MElement *ele)
{
int ndofs=0;
@@ -334,10 +365,12 @@ class xFemFunctionSpace : public FunctionSpace<T>
public:
typedef typename TensorialTraits<T>::ValType ValType;
typedef typename TensorialTraits<T>::GradType GradType;
+ typedef typename TensorialTraits<T>::HessType HessType;
protected:
FunctionSpace<T>* _spacebase;
//Function<double>* _enrichment;
public:
+ virtual void hessfuvw(MElement *ele, double u, double v, double w,std::vector<HessType> &hess);
xFemFunctionSpace(FunctionSpace<T>* spacebase) : _spacebase(spacebase) {};
virtual void f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals){};
virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads){};
@@ -354,13 +387,16 @@ class FilteredFunctionSpace : public FunctionSpace<T>
typedef typename TensorialTraits<T>::ValType ValType;
typedef typename TensorialTraits<T>::GradType GradType;
-
+ typedef typename TensorialTraits<T>::HessType HessType;
protected:
FunctionSpace<T>* _spacebase;
+
F *_filter;
public:
+
+ virtual void hessfuvw(MElement *ele, double u, double v, double w,std::vector<HessType> &hess);
FilteredFunctionSpace<T,F>(FunctionSpace<T>* spacebase,F * filter) : _spacebase(spacebase),_filter(filter)
{};
virtual void f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals){};;
diff --git a/Solver/groupOfElements.h b/Solver/groupOfElements.h
index 535861f0aa0dc6e01a1d89650bf6fd448605e76a..6a24ef7764c1dcff1b37d409df215fbb4e74ae55 100644
--- a/Solver/groupOfElements.h
+++ b/Solver/groupOfElements.h
@@ -30,14 +30,14 @@ public:
groupOfElements (GFace*);
- void addPhysical(int dim, int physical) {
+ virtual void addPhysical(int dim, int physical) {
elementFilterTrivial filter;
addPhysical (dim, physical, filter);
}
- void addElementary(GEntity *ge, const elementFilter &f);
+ virtual void addElementary(GEntity *ge, const elementFilter &f);
- void addPhysical(int dim, int physical, const elementFilter &);
+ virtual void addPhysical(int dim, int physical, const elementFilter &);
vertexContainer::const_iterator vbegin() const {
return _vertices.begin();
diff --git a/Solver/solverAlgorithms.h b/Solver/solverAlgorithms.h
index d605063fad2b357a0bbd9e7ea847a3d60e6c1464..19f21393c616a547883b23d82f71063e1507c8f0 100644
--- a/Solver/solverAlgorithms.h
+++ b/Solver/solverAlgorithms.h
@@ -1,7 +1,7 @@
//
// C++ Interface: solverAlgorithms
//
-// Description:
+// Description:
//
//
// Author: <Eric Bechet>, (C) 2009
@@ -145,7 +145,7 @@ template<class Assembler> void FixNodalDofs(FunctionSpaceBase &space,MElement *e
space.getKeys(e,R);
tabV.reserve(nv);
for (int i=0;i<nv;++i) tabV.push_back(e->getVertex(i));
-
+
for (std::vector<Dof>::iterator itd=R.begin();itd!=R.end();++itd)
{
Dof key=*itd;
diff --git a/Solver/solverField.h b/Solver/solverField.h
index 1b8e2e94841ae8e13416908359a8852bd6970b11..0160ff55ebd6b0a51fa40f8e2209eb3d07ee1122 100644
--- a/Solver/solverField.h
+++ b/Solver/solverField.h
@@ -1,7 +1,7 @@
//
// C++ Interface: field
//
-// Description:
+// Description:
//
//
// Author: <Eric Bechet>, (C) 2009
@@ -26,6 +26,7 @@ class SolverField : public FunctionSpace<T> // being able to use it instead of a
public:
typedef typename TensorialTraits<T>::ValType ValType;
typedef typename TensorialTraits<T>::GradType GradType;
+ typedef typename TensorialTraits<T>::HessType HessType;
private:
dofManager<double> *dm;
FunctionSpace<T> *fs;
@@ -71,12 +72,33 @@ class SolverField : public FunctionSpace<T> // being able to use it instead of a
grad+= SFGrads[i] * DMVals[i];
}
+ // A quoi sert cette fonction ?? (Evalue le hessien au noeuds
+/* virtual void hessf(MElement *ele, double u, double v, double w,HessType &hess)
+ {
+ // Pas besoin des dof etc
+// std::vector<Dof> D;
+ std::vector<HessType> SFHess;
+// std::vector<double> DMVals;
+// fs->getKeys(ele,D);
+// dm->getDofValue(D,DMVals);
+ fs->hessf(ele,u,v,w,SFHess);
+
+// hess=HessType;
+// for (int i=0;i<D.size();++i)
+// hess+= SFHess[i] * DMVals[i];
+ }*/
+
virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
{
GradType grad;
gradf(ele,u,v,w,grad);
grads.push_back(grad);
}
+ virtual void hessfuvw(MElement *ele, double u, double v, double w,std::vector<HessType> &hess)
+ {
+ //HessType hes;
+ fs->hessfuvw(ele,u,v,w,hess);
+ }
};
/*
@@ -87,7 +109,7 @@ class Formulation
std::vector<FunctionSpace<SVector3>* > vectorfs;
std::vector<groupOfElements* > groups;
std::vector<std::pair<MElement*,std::vector<groupOfElements&> > > links;
- dofManager<double> *dm; //
+ dofManager<double> *dm; //
};
*/
diff --git a/Solver/terms.h b/Solver/terms.h
index 4c9d9fd1a554c1c37eff34b7c37ee8e2d317aa19..36cc38f7cca5867fd7685d0396a7941fbe18f2b9 100644
--- a/Solver/terms.h
+++ b/Solver/terms.h
@@ -1,7 +1,7 @@
//
// C++ Interface: terms
//
-// Description:
+// Description:
//
//
// Author: <Eric Bechet>, (C) 2009
@@ -21,12 +21,17 @@
#include "functionSpace.h"
#include "groupOfElements.h"
#include "materialLaw.h"
-
+#include "../contrib/cm3/MInterfaceElement.h"
class BilinearTermBase
{
public :
virtual ~BilinearTermBase() {}
- virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m) =0;
+ virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m) = 0 ;
+ 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
@@ -57,7 +62,7 @@ template<class T1> class LinearTerm : public LinearTermBase
class ScalarTermBase
{
- public :
+ public :
virtual ~ScalarTermBase() {}
virtual void get(MElement *ele,int npts,IntPt *GP,double &val) =0;
};
@@ -109,9 +114,9 @@ class ScalarTermConstant : public ScalarTerm
-template<class T1,class T2> class LaplaceTerm : public BilinearTerm<T1,T2>
+template<class T1,class T2> class LaplaceTerm : public BilinearTerm<T1,T2>
{
- public :
+ public :
LaplaceTerm(FunctionSpace<T1>& space1_,FunctionSpace<T2>& space2_) : BilinearTerm<T1,T2>(space1_,space2_)
{}
virtual ~LaplaceTerm() {}
@@ -129,7 +134,7 @@ template<class T1,class T2> class LaplaceTerm : public BilinearTerm<T1,T2>
template<class T1> class LaplaceTerm<T1,T1> : public BilinearTerm<T1,T1> // symmetric
{
- public :
+ public :
LaplaceTerm(FunctionSpace<T1>& space1_) : BilinearTerm<T1,T1>(space1_,space1_)
{}
virtual ~LaplaceTerm() {}
@@ -163,9 +168,9 @@ template<class T1> class LaplaceTerm<T1,T1> : public BilinearTerm<T1,T1> // symm
-class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3>
+class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3>
{
- protected :
+ protected :
double E,nu;
bool sym;
fullMatrix<double> H;/* =
@@ -176,7 +181,7 @@ class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3>
{ 0, 0, 0, 0, C44, 0},
{ 0, 0, 0, 0, 0, C44} };*/
- public :
+ public :
IsotropicElasticTerm(FunctionSpace<SVector3>& space1_,FunctionSpace<SVector3>& space2_,double E_,double nu_) : BilinearTerm<SVector3,SVector3>(space1_,space2_),E(E_),nu(nu_),H(6,6)
{
double FACT = E / (1 + nu);
@@ -209,8 +214,10 @@ class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3>
fullMatrix<double> B(6, nbFF);
fullMatrix<double> BTH(nbFF, 6);
fullMatrix<double> BT(nbFF, 6);
+ printf("npts : %d\n",npts);
m.resize(nbFF, nbFF);
m.setAll(0.);
+ std::cout << m.size1() << " " << m.size2() << std::endl;
for (int i = 0; i < npts; i++)
{
const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2];
@@ -242,6 +249,7 @@ class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3>
fullMatrix<double> BT(nbFF1, 6);
m.resize(nbFF1, nbFF2);
m.setAll(0.);
+ // Sum on Gauss Points i
for (int i = 0; i < npts; i++)
{
const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2];
@@ -270,13 +278,13 @@ class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3>
}
BTH.setAll(0.);
BTH.gemm(BT, H);
+ // gemm add the product to m so there is a sum on gauss' points here
m.gemm(BTH, B, weight * detJ, 1.);
}
}
}
}; // class
-
inline double dot(const double &a, const double &b)
{ return a*b; }
@@ -284,7 +292,7 @@ inline double dot(const double &a, const double &b)
template<class T1> class LoadTerm : public LinearTerm<T1>
{
simpleFunction<typename TensorialTraits<T1>::ValType> &Load;
- public :
+ public :
LoadTerm(FunctionSpace<T1>& space1_,simpleFunction<typename TensorialTraits<T1>::ValType> &Load_) :LinearTerm<T1>(space1_),Load(Load_) {}
virtual ~LoadTerm() {}
diff --git a/contrib/cm3/AssembleInterface.h b/contrib/cm3/AssembleInterface.h
new file mode 100644
index 0000000000000000000000000000000000000000..be1e55d9e8378d7b7421aca70c29ba9227632fb2
--- /dev/null
+++ b/contrib/cm3/AssembleInterface.h
@@ -0,0 +1,54 @@
+//
+// 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
new file mode 100644
index 0000000000000000000000000000000000000000..5048a2b8577f7362b38fe7b0e62ada7c056e15ad
--- /dev/null
+++ b/contrib/cm3/C0DgPlateTerms.h
@@ -0,0 +1,650 @@
+//
+// 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
new file mode 100644
index 0000000000000000000000000000000000000000..8913c4390af9799b67f5ba10f872c2ad443f5615
--- /dev/null
+++ b/contrib/cm3/CMakeLists.txt
@@ -0,0 +1,8 @@
+# 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/CMakeLists.txt~ b/contrib/cm3/CMakeLists.txt~
new file mode 100644
index 0000000000000000000000000000000000000000..ecca2e400fe1669cf2a821dbe7a2c285985de57a
--- /dev/null
+++ b/contrib/cm3/CMakeLists.txt~
@@ -0,0 +1,8 @@
+# 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/GModelDg.cpp contrib/cm3/C0DgPlateTerms.h ${GMSH_SRC})
+target_link_libraries(dgsolver ${LINK_LIBRARIES})
+
diff --git a/contrib/cm3/DGterms.h~ b/contrib/cm3/DGterms.h~
new file mode 100644
index 0000000000000000000000000000000000000000..de333116512aef2fe322a3b97d9687cfd84be4c9
--- /dev/null
+++ b/contrib/cm3/DGterms.h~
@@ -0,0 +1,331 @@
+//
+// C++ Interface: terms
+//
+// Description:
+//
+//
+// Author: <Eric Bechet>, (C) 2009
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+
+#ifndef _TERMS_H_
+#define _TERMS_H_
+
+#include "SVector3.h"
+#include <vector>
+#include <iterator>
+#include "Numeric.h"
+#include "functionSpace.h"
+#include "groupOfElements.h"
+#include "materialLaw.h"
+
+class BilinearTermBase
+{
+ public :
+ virtual ~BilinearTermBase() {}
+ virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m) =0;
+};
+
+template<class T1,class T2> class BilinearTerm : public BilinearTermBase
+{
+ protected :
+ FunctionSpace<T1>& space1;
+ FunctionSpace<T2>& space2;
+ public :
+ BilinearTerm(FunctionSpace<T1>& space1_,FunctionSpace<T1>& space2_) : space1(space1_),space2(space2_) {}
+ virtual ~BilinearTerm() {}
+};
+
+class LinearTermBase
+{
+ public:
+ virtual ~LinearTermBase() {}
+ virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &v) =0;
+ virtual void get(MVertex *ver,fullVector<double> &m) =0;
+};
+
+template<class T1> class LinearTerm : public LinearTermBase
+{
+ protected :
+ FunctionSpace<T1>& space1;
+ public :
+ LinearTerm(FunctionSpace<T1>& space1_) : space1(space1_) {}
+ virtual ~LinearTerm() {}
+};
+
+class ScalarTermBase
+{
+ public :
+ virtual ~ScalarTermBase() {}
+ virtual void get(MElement *ele,int npts,IntPt *GP,double &val) =0;
+};
+
+class ScalarTerm : public ScalarTermBase
+{
+ public :
+ virtual ~ScalarTerm() {}
+};
+
+
+template<class T1,class T2> class BilinearTermToScalarTerm : public ScalarTerm
+{
+ BilinearTerm<T1,T2> &bilterm;
+ public :
+ BilinearTermToScalarTerm(BilinearTerm<T1,T2> &bilterm_): bilterm(bilterm_){}
+ virtual ~BilinearTermToScalarTerm() {}
+ virtual void get(MElement *ele,int npts,IntPt *GP,double &val)
+ {
+ fullMatrix<double> localMatrix;
+ bilterm.get(ele,npts,GP,localMatrix);
+ val=localMatrix(0,0);
+ }
+};
+
+
+class ScalarTermConstant : public ScalarTerm
+{
+ double val;
+ public :
+ ScalarTermConstant(double val_=1.0): val(val_) {}
+ virtual ~ScalarTermConstant() {}
+ virtual void get(MElement *ele,int npts,IntPt *GP,double &val)
+ {
+ double jac[3][3];
+ val=0;
+ for (int i = 0; i < npts; 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 = ele->getJacobian(u, v, w, jac);
+ val+=weight*detJ;
+ }
+ }
+ virtual void get(MVertex *ver,double &val)
+ {
+ val=1;
+ }
+};
+
+
+
+template<class T1,class T2> class LaplaceTerm : public BilinearTerm<T1,T2>
+{
+ public :
+ LaplaceTerm(FunctionSpace<T1>& space1_,FunctionSpace<T2>& space2_) : BilinearTerm<T1,T2>(space1_,space2_)
+ {}
+ virtual ~LaplaceTerm() {}
+ virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m)
+ {
+ Msg::Error("LaplaceTerm<S1,S2> w/ S1!=S2 not implemented");
+ }
+ virtual void get(MVertex *ver,fullMatrix<double> &m)
+ {
+ Msg::Error("LaplaceTerm<S1,S2> w/ S1!=S2 not implemented");
+ }
+}; // class
+
+
+
+template<class T1> class LaplaceTerm<T1,T1> : public BilinearTerm<T1,T1> // symmetric
+{
+ public :
+ LaplaceTerm(FunctionSpace<T1>& space1_) : BilinearTerm<T1,T1>(space1_,space1_)
+ {}
+ virtual ~LaplaceTerm() {}
+ virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m)
+ {
+ int nbFF = BilinearTerm<T1,T1>::space1.getNumKeys(ele);
+ double jac[3][3];
+ m.resize(nbFF, nbFF);
+ m.setAll(0.);
+ for (int i = 0; i < npts; 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 = ele->getJacobian(u, v, w, jac);
+ std::vector<typename TensorialTraits<T1>::GradType> Grads;
+ BilinearTerm<T1,T1>::space1.gradf(ele,u, v, w, Grads);
+ for (int j = 0; j < nbFF; j++)
+ {
+ for (int k = j; k < nbFF; k++)
+ {
+ double contrib=weight * detJ * dot(Grads[j],Grads[k]);
+ m(j,k)+=contrib;
+ if (j!=k) m(k,j)+=contrib;
+ }
+ }
+ }
+// m.print("");
+// exit(0);
+ }
+}; // class
+
+
+
+
+class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3>
+{
+ protected :
+ double E,nu;
+ bool sym;
+ fullMatrix<double> H;/* =
+ { {C11, C12, C12, 0, 0, 0},
+ {C12, C11, C12, 0, 0, 0},
+ {C12, C12, C11, 0, 0, 0},
+ { 0, 0, 0, C44, 0, 0},
+ { 0, 0, 0, 0, C44, 0},
+ { 0, 0, 0, 0, 0, C44} };*/
+
+ public :
+ IsotropicElasticTerm(FunctionSpace<SVector3>& space1_,FunctionSpace<SVector3>& space2_,double E_,double nu_) : BilinearTerm<SVector3,SVector3>(space1_,space2_),E(E_),nu(nu_),H(6,6)
+ {
+ double FACT = E / (1 + nu);
+ double C11 = FACT * (1 - nu) / (1 - 2 * nu);
+ double C12 = FACT * nu / (1 - 2 * nu);
+ double C44 = (C11 - C12) / 2;
+ H.scale(0.);
+ for (int i=0;i<3;++i) {H(i,i)=C11;H(i+3,i+3)=C44;}
+ H(1,0)=H(0,1)=H(2,0)=H(0,2)=H(1,2)=H(2,1)=C12;
+ sym=(&space1_==&space2_);
+ }
+ IsotropicElasticTerm(FunctionSpace<SVector3>& space1_,double E_,double nu_) : BilinearTerm<SVector3,SVector3>(space1_,space1_),E(E_),nu(nu_),H(6,6)
+ {
+ double FACT = E / (1 + nu);
+ double C11 = FACT * (1 - nu) / (1 - 2 * nu);
+ double C12 = FACT * nu / (1 - 2 * nu);
+ double C44 = (C11 - C12) / 2;
+ H.scale(0.);
+ for (int i=0;i<3;++i) {H(i,i)=C11;H(i+3,i+3)=C44;}
+ H(1,0)=H(0,1)=H(2,0)=H(0,2)=H(1,2)=H(2,1)=C12;
+ sym=true;
+ }
+ virtual ~IsotropicElasticTerm() {}
+ virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m)
+ {
+ if (sym)
+ {
+ int nbFF = BilinearTerm<SVector3,SVector3>::space1.getNumKeys(ele);
+ double jac[3][3];
+ fullMatrix<double> B(6, nbFF);
+ fullMatrix<double> BTH(nbFF, 6);
+ fullMatrix<double> BT(nbFF, 6);
+ m.resize(nbFF, nbFF);
+ m.setAll(0.);
+ for (int i = 0; i < npts; 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 = ele->getJacobian(u, v, w, jac);
+ std::vector<TensorialTraits<SVector3>::GradType> Grads;
+ BilinearTerm<SVector3,SVector3>::space1.gradf(ele,u, v, w, Grads); // a optimiser ??
+ for (int j = 0; j < nbFF; j++)
+ {
+ BT(j, 0) = B(0, j) = Grads[j](0,0);
+ BT(j, 1) = B(1, j) = Grads[j](1,1);
+ BT(j, 2) = B(2, j) = Grads[j](2,2);
+ BT(j, 3) = B(3, j) = Grads[j](0,1)+Grads[j](1,0);
+ BT(j, 4) = B(4, j) = Grads[j](1,2)+Grads[j](2,1);
+ BT(j, 5) = B(5, j) = Grads[j](0,2)+Grads[j](2,0);
+ }
+ BTH.setAll(0.);
+ BTH.gemm(BT, H);
+ m.gemm(BTH, B, weight * detJ, 1.);
+ }
+ }
+ else
+ {
+ int nbFF1 = BilinearTerm<SVector3,SVector3>::space1.getNumKeys(ele);
+ int nbFF2 = BilinearTerm<SVector3,SVector3>::space2.getNumKeys(ele);
+ double jac[3][3];
+ fullMatrix<double> B(6, nbFF2);
+ fullMatrix<double> BTH(nbFF2, 6);
+ fullMatrix<double> BT(nbFF1, 6);
+ m.resize(nbFF1, nbFF2);
+ m.setAll(0.);
+ for (int i = 0; i < npts; 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 = ele->getJacobian(u, v, w, jac);
+ std::vector<TensorialTraits<SVector3>::GradType> Grads;// tableau de matrices...
+ std::vector<TensorialTraits<SVector3>::GradType> GradsT;// tableau de matrices...
+ BilinearTerm<SVector3,SVector3>::space1.gradf(ele,u, v, w, Grads);
+ BilinearTerm<SVector3,SVector3>::space2.gradf(ele,u, v, w, GradsT);
+ for (int j = 0; j < nbFF1; j++)
+ {
+ BT(j, 0) = Grads[j](0,0);
+ BT(j, 1) = Grads[j](1,1);
+ BT(j, 2) = Grads[j](2,2);
+ BT(j, 3) = Grads[j](0,1)+Grads[j](1,0);
+ BT(j, 4) = Grads[j](1,2)+Grads[j](2,1);
+ BT(j, 5) = Grads[j](0,2)+Grads[j](2,0);
+ }
+ for (int j = 0; j < nbFF2; j++)
+ {
+ B(0, j) = GradsT[j](0,0);
+ B(1, j) = GradsT[j](1,1);
+ B(2, j) = GradsT[j](2,2);
+ B(3, j) = GradsT[j](0,1)+GradsT[j](1,0);
+ B(4, j) = GradsT[j](1,2)+GradsT[j](2,1);
+ B(5, j) = GradsT[j](0,2)+GradsT[j](2,0);
+ }
+ BTH.setAll(0.);
+ BTH.gemm(BT, H);
+ m.gemm(BTH, B, weight * detJ, 1.);
+ }
+ }
+ }
+}; // class
+
+
+inline double dot(const double &a, const double &b)
+{ return a*b; }
+
+
+template<class T1> class LoadTerm : public LinearTerm<T1>
+{
+ simpleFunction<typename TensorialTraits<T1>::ValType> &Load;
+ public :
+ LoadTerm(FunctionSpace<T1>& space1_,simpleFunction<typename TensorialTraits<T1>::ValType> &Load_) :LinearTerm<T1>(space1_),Load(Load_) {}
+ virtual ~LoadTerm() {}
+
+ virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &m)
+ {
+ double nbFF=LinearTerm<T1>::space1.getNumKeys(ele);
+ double jac[3][3];
+ m.resize(nbFF);
+ m.scale(0.);
+ for (int i = 0; i < npts; 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 = ele->getJacobian(u, v, w, jac);
+ std::vector<typename TensorialTraits<T1>::ValType> Vals;
+ LinearTerm<T1>::space1.f(ele,u, v, w, Vals);
+ SPoint3 p;
+ ele->pnt(u, v, w, p);
+ typename TensorialTraits<T1>::ValType load=Load(p.x(),p.y(),p.z());
+ for (int j = 0; j < nbFF ; ++j)
+ {
+ m(j)+=dot(Vals[j],load)*weight*detJ;
+ }
+ }
+ }
+
+ virtual void get(MVertex *ver,fullVector<double> &m)
+ {
+ double nbFF=LinearTerm<T1>::space1.getNumKeys(ver);
+ double jac[3][3];
+ m.resize(nbFF);
+ std::vector<typename TensorialTraits<T1>::ValType> Vals;
+ LinearTerm<T1>::space1.f(ver, Vals);
+ typename TensorialTraits<T1>::ValType load=Load(ver->x(),ver->y(),ver->z());
+ for (int j = 0; j < nbFF ; ++j)
+ {
+ m(j)=dot(Vals[j],load);
+ }
+ }
+};
+
+
+
+
+#endif// _TERMS_H_
diff --git a/contrib/cm3/DgC0PlateElementaryTerms.h b/contrib/cm3/DgC0PlateElementaryTerms.h
new file mode 100644
index 0000000000000000000000000000000000000000..b79a215e99495693dd55124d59ef71474ff8a462
--- /dev/null
+++ b/contrib/cm3/DgC0PlateElementaryTerms.h
@@ -0,0 +1,387 @@
+//
+// 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/DgC0PlateElementaryTerms1ddl.h~ b/contrib/cm3/DgC0PlateElementaryTerms1ddl.h~
new file mode 100644
index 0000000000000000000000000000000000000000..1fa6c24644a21068f628fa193904356573e699b6
--- /dev/null
+++ b/contrib/cm3/DgC0PlateElementaryTerms1ddl.h~
@@ -0,0 +1,69 @@
+// Compute the component (j,k) of the elementary stiffness matrix
+inline double BulkC0PlateDGStiffnessBendingTerms(TensorialTraits<double>::HessType &hessj, TensorialTraits<double>::HessType &hessk, LinearElasticShellHookeTensor *H, const LocalBasis *lb){
+ double val = 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++)
+ val += hessj(alpha,beta)*hessk(gamma,delta)*H->get(alpha,beta,gamma,delta);
+ return lb->gett0(2)*lb->gett0(2)*val;
+}
+
+inline double consAndCompC0PlateStiffnessTerms(LinearElasticShellHookeTensor *Hhat,const fullMatrix<double> &Bhat,const SVector3 &dt, const LocalBasis *lb){
+ double val=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++)
+ val += Hhat->get(alpha,beta,gamma,delta)*Bhat(gamma,delta)*dot(dt,lb->getphi0(alpha))*(-lb->getphi0(1,beta));
+ return 0.5*val;
+}
+
+inline double stabilityC0PlateStiffnessTerms(LinearElasticShellHookeTensor *Hhat, const SVector3 &dta, const SVector3 &dtb, const LocalBasis *lb){
+ double val=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++)
+ val += dot(dta,lb->getphi0(gamma))*dot(dtb,lb->getphi0(alpha))*Hhat->get(alpha,beta,gamma,delta)*(-lb->getphi0(1,delta))*(-lb->getphi0(1,beta));
+
+ return val;
+}
+
+inline double BulkC0PlateDGForceTerms(const TensorialTraits<double>::HessType &hessj,const std::vector<TensorialTraits<double>::HessType> &Hess,const LinearElasticShellHookeTensor *H,const LocalBasis *lb,const fullMatrix<double> &disp){
+ const int n = Hess.size();
+ double sum,val;
+ val=0.;
+ for(int a=0;a<2;a++)
+ for(int b=0;b<2;b++){
+ sum=0.;
+ for(int j=0;j<n;j++)
+ sum+=-lb->gett0(2)*Hess[j](a,b)*disp(j,0);
+ for(int c=0;c<2;c++)
+ for(int d=0;d<2;d++)
+ val += -lb->gett0(2)*hessj(c,d)*H->get(a,b,c,d)*sum;
+ }
+ return val;
+}
+
+inline double consC0PlateForceTerms(const LinearElasticShellHookeTensor *Hhat, const fullMatrix<double> &Bhat, const std::vector<SVector3> &Dt_m, const std::vector<SVector3> &Dt_p,const LocalBasis *lb, const fullMatrix<double> &disp){
+ const int n_m = Dt_m.size();
+ const int n_p = Dt_p.size();
+ double sum,val;
+ val=0.;
+ for(int a=0;a<2;a++){
+ sum=0.;
+ for(int j=0;j<n_m;j++)
+ sum -= dot(Dt_m[j],lb->getphi0(a))*disp(j,0);
+ for(int j=0;j<n_p;j++)
+ sum += dot(Dt_p[j],lb->getphi0(a))*disp(j+n_m,0);
+
+ for(int b=0;b<2;b++)
+ for(int c=0;c<2;c++)
+ for(int d=0;d<2;d++)
+ val+= Hhat->get(a,b,c,d)*Bhat(c,d)*sum*(-lb->getphi0(1,b));
+ }
+ return 0.5*val;
+}
+
+
diff --git a/contrib/cm3/DgC0PlateSolver.cpp b/contrib/cm3/DgC0PlateSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b4bbd36a29a1475ea48286c00e50457fd51848b0
--- /dev/null
+++ b/contrib/cm3/DgC0PlateSolver.cpp
@@ -0,0 +1,386 @@
+// 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
new file mode 100644
index 0000000000000000000000000000000000000000..0d3689426b873ad5e4b447bb87e0e72308bb63c4
--- /dev/null
+++ b/contrib/cm3/DgC0PlateSolver.h
@@ -0,0 +1,96 @@
+// 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
new file mode 100644
index 0000000000000000000000000000000000000000..cad85d7349b7b442d9a41098dd32aa8eb66c5bad
--- /dev/null
+++ b/contrib/cm3/GModelWithInterface.cpp
@@ -0,0 +1,664 @@
+//
+// 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
new file mode 100644
index 0000000000000000000000000000000000000000..ad07ff51711a8533a6307616a7b21994b7c294d3
--- /dev/null
+++ b/contrib/cm3/GModelWithInterface.h
@@ -0,0 +1,59 @@
+//
+// 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
new file mode 100644
index 0000000000000000000000000000000000000000..e00a2dbbae4f5bff50f8b21f0b65c560bda0e6dc
--- /dev/null
+++ b/contrib/cm3/LinearElasticShellHookeTensor.h
@@ -0,0 +1,75 @@
+//
+// 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
new file mode 100644
index 0000000000000000000000000000000000000000..e39f24e93ab600a1acab5e033e49778f33dcf528
--- /dev/null
+++ b/contrib/cm3/LocalBasis.cpp
@@ -0,0 +1,13 @@
+//
+// 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
new file mode 100644
index 0000000000000000000000000000000000000000..04ede3ff17640f149bcb24c031fcb92ac271ebcf
--- /dev/null
+++ b/contrib/cm3/LocalBasis.h
@@ -0,0 +1,174 @@
+//
+// 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
new file mode 100644
index 0000000000000000000000000000000000000000..3c5e97520ed156c5e3030a8a57d4f9f5c9e8b213
--- /dev/null
+++ b/contrib/cm3/MInterfaceElement.cpp
@@ -0,0 +1,13 @@
+//
+// 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
new file mode 100644
index 0000000000000000000000000000000000000000..a2c20e23076343079c4bcf477758ea8a747f5b65
--- /dev/null
+++ b/contrib/cm3/MInterfaceElement.h
@@ -0,0 +1,167 @@
+//
+// 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(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(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
new file mode 100644
index 0000000000000000000000000000000000000000..af5627e118eb36009df212b52fe54a8fa339a98c
--- /dev/null
+++ b/contrib/cm3/Tensor4dim2.h
@@ -0,0 +1,116 @@
+//
+// 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
new file mode 100644
index 0000000000000000000000000000000000000000..7f8f0382a038e0ee7edabee4e05cadfbab00f561
--- /dev/null
+++ b/contrib/cm3/mainDG.cpp
@@ -0,0 +1,39 @@
+#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();
+
+}
diff --git a/contrib/cm3/mainElasticity.cpp~ b/contrib/cm3/mainElasticity.cpp~
new file mode 100644
index 0000000000000000000000000000000000000000..3ac287aedcef43dfb36e8885415709947394ebd9
--- /dev/null
+++ b/contrib/cm3/mainElasticity.cpp~
@@ -0,0 +1,108 @@
+#include "Gmsh.h"
+#include "elasticitySolver.h"
+#include "PView.h"
+#include "PViewData.h"
+#include "highlevel.h"/usr/bin/gmake -f "/home/gauthier/distgmsh/trunk/Makefile" elastic -o ../../bin/elastic
+#include "groupOfElements.h"
+#include <iterator>
+
+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
+ elasticitySolver 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");
+ pv->getData()->writeMSH("energ.msh", false);
+ delete pv;
+
+
+
+ // stop gmsh
+ GmshFinalize();
+
+}
+
+
+
+/*
+ groupOfElements *g = new groupOfElements (2, 7);
+
+ MElement *e=*(g->begin());
+ std::cout << e->getNumPrimaryVertices() << "vertices" << std::endl;
+ const double uvw[3]={0.,0.,0.};
+ std::vector<Dof> dofs;
+ std::vector<double> vals;
+ std::vector<SVector3> grads;
+ std::vector<SVector3> vals2;
+ std::vector<STensor3> grads2;
+
+ std::ostream_iterator< double > output( std::cout, " " );
+
+ ScalarLagrangeFunctionSpace L(100);
+ std::cout << L.getNumKeys(e) << "fonctions de formes L" << std::endl;
+ L.getKeys(e,dofs);
+ for (int i=0;i<dofs.size();++i) std::cout << "entity: " << dofs[i].getEntity() << " id: " << dofs[i].getType() << std::endl;
+ dofs.clear();
+ L.f(e,0.1,0.1,0,vals);
+ L.gradf(e,0.1,0.1,0,grads);
+ std::copy(vals.begin(),vals.end(),output); std::cout << std::endl;
+ for (std::vector<SVector3>::iterator it=grads.begin();it!=grads.end();++it) { std::cout << (*it)[0]<< " " << (*it)[1] <<" " << (*it)[2] << std::endl; }
+
+ VectorLagrangeFunctionSpace L1(100,VectorLagrangeFunctionSpace::VECTOR_X);
+ VectorLagrangeFunctionSpace L2(100,VectorLagrangeFunctionSpace::VECTOR_Y);
+ std::cout << L2.getNumKeys(e) << "fonctions de formes L2" << std::endl;
+ L2.f(e,0.1,0.1,0,vals2);
+ L2.gradf(e,0.1,0.1,0,grads2);
+ for (std::vector<SVector3>::iterator it=vals2.begin();it!=vals2.end();++it) { std::cout << (*it)[0]<< " " << (*it)[1] <<" " << (*it)[2] << std::endl; }
+ for (std::vector<STensor3>::iterator it=grads2.begin();it!=grads2.end();++it) { (*it).print(""); }
+
+ VectorLagrangeFunctionSpace L3(100,VectorLagrangeFunctionSpace::VECTOR_Z);
+
+ VectorLagrangeFunctionSpace P123(100);
+ std::cout << P123.getNumKeys(e) << "fonctions de formes P123" << std::endl;
+ P123.getKeys(e,dofs);
+ std::cout << dofs.size() << std::endl;
+ for (int i=0;i<dofs.size();++i) std::cout << "entity: " << dofs[i].getEntity() << " id: " << dofs[i].getType() << std::endl;
+
+ vals2.clear();
+ grads2.clear();
+ P123.f(e,0.1,0.1,0,vals2);
+ P123.gradf(e,0.1,0.1,0,grads2);
+ for (std::vector<SVector3>::iterator it=vals2.begin();it!=vals2.end();++it) { std::cout << (*it)[0]<< " " << (*it)[1] <<" " << (*it)[2] << std::endl; }
+ for (std::vector<STensor3>::iterator it=grads2.begin();it!=grads2.end();++it) { (*it).print(""); }
+
+
+
+ FormBilinear<TermBilinearMeca,ScalarLagrangeFunctionSpace,ScalarLagrangeFunctionSpace > f(L,L);
+ f.func();
+ f.Accumulate(e,uvw);
+
+
+ FormBilinear<TermBilinearMecaNL,ScalarLagrangeFunctionSpace,ScalarLagrangeFunctionSpace > fnl(L,L);
+ fnl.func();
+
+*/
+
+
+
+
+
diff --git a/utils/api_demos/Makefile b/utils/api_demos/Makefile
index 929e1bc7a0de5c7b4668283a460053d81a947908..5089196224023e66692844abd7eb23370d9c5154 100644
--- a/utils/api_demos/Makefile
+++ b/utils/api_demos/Makefile
@@ -1,6 +1,7 @@
FLAGS = -g -O2 ${CXXFLAGS}
-INC = -I/usr/local/include -I.
-GMSH = -L/usr/local/lib -lGmsh
+INC = -I/usr/local/include -I. -I /home/gauthier/distgmsh/trunk/buildlib/lib/include
+GMSH = -L/home/gauthier/distgmsh/trunk/buildlib/lib/lib -lGmsh
+METIS= -L/usr/local/lib -lmetis
GLUT = -framework OpenGL -framework GLUT -framework Cocoa -framework ApplicationServices
#GLUT = -lGLUT -lGLU -lGL -lX11
OCC = -L/usr/local/opencascade/lib -L/Users/remacle/SOURCES/opencascade/lib -lTKSTEP -lTKSTEP209 -lTKSTEPAttr -lTKSTEPBase -lTKIGES -lTKXSBase -lTKOffset -lTKFeat -lTKFillet -lTKBool -lTKShHealing -lTKMesh -lTKHLR -lTKBO -lTKPrim -lTKTopAlgo -lTKGeomAlgo -lTKBRep -lTKGeomBase -lTKG3d -lTKG2d -lTKAdvTools -lTKMath -lTKernel