Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[ContactStructuralMechanicsApplication] Adding new debugging info/pythonic clean up #3582

Merged
merged 3 commits into from
Dec 5, 2018
Merged

[ContactStructuralMechanicsApplication] Adding new debugging info/pythonic clean up #3582

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
37e9315
More debugging info into contact
Dec 4, 2018
f1c2c13
Less verbosity
Dec 4, 2018
a4d2d1b
Merge branch 'master' into contact/more-debugging-info
Dec 5, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
123 changes: 76 additions & 47 deletions applications/ContactStructuralMechanicsApplication/python_scripts/search_base_process.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,9 +85,9 @@ def __init__(self, Model, settings):
self.database_step = 0

# Detect "End" as a tag and replace it by a large number
if(self.settings.Has("interval")):
if(self.settings["interval"][1].IsString() ):
if(self.settings["interval"][1].GetString() == "End"):
if self.settings.Has("interval"):
if self.settings["interval"][1].IsString():
if self.settings["interval"][1].GetString() == "End":
self.settings["interval"][1].SetDouble(1e30) # = default_settings["interval"][1]
else:
raise Exception("the second value of interval can be \"End\" or a number, interval currently:"+settings["interval"].PrettyPrintJsonString())
Expand All @@ -103,7 +103,7 @@ def ExecuteInitialize(self):
"""

# First we generate or identify the different model parts
if (self.computing_model_part.HasSubModelPart("Contact")):
if self.computing_model_part.HasSubModelPart("Contact"):
self.preprocess = False
# We get the submodelpart
self.search_model_part = self.computing_model_part.GetSubModelPart("Contact")
Expand All @@ -113,25 +113,25 @@ def ExecuteInitialize(self):
self.search_model_part = self.computing_model_part.CreateSubModelPart("Contact")

# In case of no "Contact" model part we create it
if (self.preprocess is True):
if self.preprocess is True:
# In case no model part is assigned we detect the skin
self.count_search_model_part = 0
for key in self.settings["search_model_part"].keys():
if (self.settings["search_model_part"][key].size() > 0):
if self.settings["search_model_part"][key].size() > 0:
self.count_search_model_part += 1
if (self.count_search_model_part == 0):
if self.count_search_model_part == 0:
self.__detect_skin(self.main_model_part)
else:
for key in self.settings["search_model_part"].keys():
if (self.settings["search_model_part"][key].size() > 0):
if self.settings["search_model_part"][key].size() > 0:
self.__generate_search_model_part_from_input_list(self.settings["search_model_part"][key], key)

# We compute NODAL_H that can be used in the search and some values computation
self.find_nodal_h = KM.FindNodalHProcess(self.computing_model_part)
self.find_nodal_h.Execute()

## We recompute the search factor and the check in function of the relative size of the mesh
if (self.settings["search_parameters"]["adapt_search"].GetBool() is True):
if self.settings["search_parameters"]["adapt_search"].GetBool() is True:
factor = CSMA.ContactUtilities.CalculateRelativeSizeMesh(self.computing_model_part)
KM.Logger.PrintWarning("SEARCH ADAPT FACTOR: ", "{:.2e}".format(factor))
search_factor = self.settings["search_parameters"]["search_factor"].GetDouble() * factor
Expand All @@ -143,7 +143,7 @@ def ExecuteInitialize(self):
self._initialize_process_info()

#If the conditions doesn't exist we create them
if (self.preprocess is False):
if self.preprocess is False:
master_slave_process = CSMA.MasterSlaveProcess(self.computing_model_part)
master_slave_process.Execute()

Expand All @@ -160,15 +160,15 @@ def ExecuteInitialize(self):
# Creating the search
self.search_search = {}
for key in self.settings["search_model_part"].keys():
if (self.settings["search_model_part"][key].size() > 0):
if self.settings["search_model_part"][key].size() > 0:
self._create_main_search(key)

# We initialize conditions
self._initialize_search_conditions()

# We perform search initializations
for key in self.settings["search_model_part"].keys():
if (self.settings["search_model_part"][key].size() > 0):
if self.settings["search_model_part"][key].size() > 0:
# We initialize the search utility
self.search_search[key].CheckContactModelParts()
self.search_search[key].CreatePointListMortar()
Expand All @@ -189,16 +189,16 @@ def ExecuteInitializeSolutionStep(self):
self -- It signifies an instance of a class.
"""

if(self._get_if_is_interval() is True):
if self._get_if_is_interval() is True:
self.database_step += 1
global_step = self.main_model_part.ProcessInfo[KM.STEP]
database_step_update = self.settings["search_parameters"]["database_step_update"].GetInt()
if (self.database_step >= database_step_update or global_step == 1):
if self.database_step >= database_step_update or global_step == 1:
# We unset the flag MARKER (used in the nodes to not deactivate it)
KM.VariableUtils().SetFlag(KM.MARKER, False, self._get_process_model_part().Nodes)
# We solve one linear step with a linear strategy if needed
for key in self.settings["search_model_part"].keys():
if (self.settings["search_model_part"][key].size() > 0):
if self.settings["search_model_part"][key].size() > 0:
# Clear current pairs
self.search_search[key].ClearMortarConditions()
# Update database
Expand All @@ -209,7 +209,7 @@ def ExecuteInitializeSolutionStep(self):
KM.VariableUtils().SetFlag(KM.MARKER, False, self._get_process_model_part().Nodes)

# Debug
if (self.settings["search_parameters"]["debug_mode"].GetBool() is True):
if self.settings["search_parameters"]["debug_mode"].GetBool() is True:
self._debug_output(global_step, "", self._get_problem_name())
# We compute the total integrated area, for debugging
self.__get_integration_area()
Expand All @@ -224,7 +224,36 @@ def ExecuteFinalizeSolutionStep(self):
Keyword arguments:
self -- It signifies an instance of a class.
"""
pass

# Debug we compute if the total load corresponds with the total contact force and the reactions
if self.settings["search_parameters"]["debug_mode"].GetBool() is True:
total_load = KM.Vector(3)
total_load[0] = 0.0
total_load[1] = 0.0
total_load[2] = 0.0
total_reaction = KM.Vector(3)
total_reaction[0] = 0.0
total_reaction[1] = 0.0
total_reaction[2] = 0.0
total_contact_force = 0

# Computing total load applied (I will consider only surface loads for now)
for cond in self.computing_model_part.Conditions:
geom = cond.GetGeometry()
if cond.Has(SMA.LINE_LOAD):
total_load += geom.Length() * cond.GetValue(SMA.LINE_LOAD)
if cond.Has(SMA.SURFACE_LOAD):
total_load += geom.Area() * cond.GetValue(SMA.SURFACE_LOAD)

for node in self.computing_model_part.Nodes:
if node.Has(KM.NODAL_AREA) and node.Has(CSMA.AUGMENTED_NORMAL_CONTACT_PRESSURE):
total_contact_force += node.GetValue(KM.NODAL_AREA) * node.GetValue(CSMA.AUGMENTED_NORMAL_CONTACT_PRESSURE)

total_reaction += node.GetSolutionStepValue(KM.REACTION)

KM.Logger.PrintWarning("TOTAL LOAD: ", "X: {:.2e}".format(total_load[0]), "\t Y: {:.2e}".format(total_load[1]), "\tZ: {:.2e}".format(total_load[2]))
KM.Logger.PrintWarning("TOTAL REACTION: ", "X: {:.2e}".format(total_reaction[0]), "\t Y: {:.2e}".format(total_reaction[1]), "\tZ: {:.2e}".format(total_reaction[2]))
KM.Logger.PrintWarning("TOTAL CONTACT FORCE: ", "{:.2e}".format(total_contact_force))

def ExecuteBeforeOutputStep(self):
""" This method is executed right before the ouput process computation
Expand All @@ -241,12 +270,12 @@ def ExecuteAfterOutputStep(self):
self -- It signifies an instance of a class.
"""
global_step = self.main_model_part.ProcessInfo[KM.STEP]
if(self._get_if_is_interval() is True):
if self._get_if_is_interval() is True:
modified = self.main_model_part.Is(KM.MODIFIED)
database_step_update = self.settings["search_parameters"]["database_step_update"].GetInt()
if (modified is False and (self.database_step >= database_step_update or global_step == 1)):
if modified is False and (self.database_step >= database_step_update or global_step == 1):
for key in self.settings["search_model_part"].keys():
if (self.settings["search_model_part"][key].size() > 0):
if self.settings["search_model_part"][key].size() > 0:
self.search_search[key].ClearMortarConditions()
self.database_step = 0

Expand Down Expand Up @@ -293,11 +322,11 @@ def _assign_master_flags(self, partial_model_part):
self -- It signifies an instance of a class.
"""

if (self.predefined_master_slave is True):
if self.predefined_master_slave is True:
KM.VariableUtils().SetFlag(KM.SLAVE, False, partial_model_part.Nodes)
KM.VariableUtils().SetFlag(KM.MASTER, True, partial_model_part.Nodes)

if (len(partial_model_part.Conditions) > 0):
if len(partial_model_part.Conditions) > 0:
KM.VariableUtils().SetFlag(KM.SLAVE, False, partial_model_part.Conditions)
KM.VariableUtils().SetFlag(KM.MASTER, True, partial_model_part.Conditions)

Expand All @@ -308,7 +337,7 @@ def _assign_slave_flags(self, key = "0"):
self -- It signifies an instance of a class.
"""

if (self.predefined_master_slave is True):
if self.predefined_master_slave is True:
if not self.settings["assume_master_slave"][key].IsArray():
raise Exception("{0} Error: Model part list is unreadable".format(self.__class__.__name__))
for i in range(0, self.settings["assume_master_slave"][key].size()):
Expand All @@ -317,7 +346,7 @@ def _assign_slave_flags(self, key = "0"):
KM.VariableUtils().SetFlag(KM.SLAVE, True, model_part_slave.Nodes)
KM.VariableUtils().SetFlag(KM.MASTER, False, model_part_slave.Nodes)

if (len(model_part_slave.Conditions) > 0):
if len(model_part_slave.Conditions) > 0:
KM.VariableUtils().SetFlag(KM.SLAVE, True, model_part_slave.Conditions)
KM.VariableUtils().SetFlag(KM.MASTER, False, model_part_slave.Conditions)

Expand All @@ -335,7 +364,7 @@ def _interface_preprocess(self, partial_model_part, key = "0"):
# It should create the conditions automatically
interface_parameters = KM.Parameters("""{"simplify_geometry": false, "contact_property_id": 0}""")
interface_parameters["contact_property_id"].SetInt(self.settings["search_property_ids"][key].GetInt())
if (self.dimension == 2):
if self.dimension == 2:
self.interface_preprocess.GenerateInterfacePart2D(partial_model_part, interface_parameters)
else:
self.interface_preprocess.GenerateInterfacePart3D(partial_model_part, interface_parameters)
Expand Down Expand Up @@ -409,17 +438,17 @@ def _create_main_search(self, key = "0"):
number_nodes, number_nodes_master = self._compute_number_nodes()

# We create the search process
if (self.dimension == 2):
if self.dimension == 2:
self.search_search[key] = CSMA.TreeContactSearch2D2N(self.computing_model_part, search_parameters)
else:
if (number_nodes == 3):
if (number_nodes_master == 3):
if number_nodes == 3:
if number_nodes_master == 3:
self.search_search[key] = CSMA.TreeContactSearch3D3N(self.computing_model_part, search_parameters)
else:
self.search_search[key] = CSMA.TreeContactSearch3D3N4N(self.computing_model_part, search_parameters)

elif (number_nodes == 4):
if (number_nodes_master == 3):
elif number_nodes == 4:
if number_nodes_master == 3:
self.search_search[key] = CSMA.TreeContactSearch3D4N3N(self.computing_model_part, search_parameters)
else:
self.search_search[key] = CSMA.TreeContactSearch3D4N(self.computing_model_part, search_parameters)
Expand Down Expand Up @@ -473,20 +502,20 @@ def _debug_output(self, label, name, problem_name = ""):
gid_io.WriteNodalResults(KM.NORMAL, self.main_model_part.Nodes, label, 0)
gid_io.WriteNodalResultsNonHistorical(KM.NODAL_AREA, self.main_model_part.Nodes, label)
gid_io.WriteNodalResults(KM.DISPLACEMENT, self.main_model_part.Nodes, label, 0)
if (self.main_model_part.Nodes[1].SolutionStepsDataHas(KM.VELOCITY_X) is True):
if self.main_model_part.Nodes[1].SolutionStepsDataHas(KM.VELOCITY_X) is True:
gid_io.WriteNodalResults(KM.VELOCITY, self.main_model_part.Nodes, label, 0)
gid_io.WriteNodalResults(KM.ACCELERATION, self.main_model_part.Nodes, label, 0)
gid_io.WriteNodalResultsNonHistorical(CSMA.NORMAL_GAP, self.main_model_part.Nodes, label)
gid_io.WriteNodalResultsNonHistorical(CSMA.AUXILIAR_COORDINATES, self.main_model_part.Nodes, label)
gid_io.WriteNodalResultsNonHistorical(CSMA.DELTA_COORDINATES, self.main_model_part.Nodes, label)

if (problem_name == "Contact"):
if problem_name == "Contact":
gid_io.WriteNodalResults(CSMA.WEIGHTED_GAP, self.main_model_part.Nodes, label, 0)
gid_io.WriteNodalResultsNonHistorical(CSMA.AUGMENTED_NORMAL_CONTACT_PRESSURE, self.main_model_part.Nodes, label)
if (self.is_frictional is True):
if self.is_frictional is True:
gid_io.WriteNodalFlags(KM.SLIP, "SLIP", self.main_model_part.Nodes, label)
gid_io.WriteNodalResultsNonHistorical(CSMA.AUGMENTED_TANGENT_CONTACT_PRESSURE, self.main_model_part.Nodes, label)
if (self.main_model_part.Nodes[1].SolutionStepsDataHas(CSMA.LAGRANGE_MULTIPLIER_CONTACT_PRESSURE) is True):
if self.main_model_part.Nodes[1].SolutionStepsDataHas(CSMA.LAGRANGE_MULTIPLIER_CONTACT_PRESSURE) is True:
gid_io.WriteNodalResults(CSMA.LAGRANGE_MULTIPLIER_CONTACT_PRESSURE, self.main_model_part.Nodes, label, 0)
else:
gid_io.WriteNodalResults(KM.VECTOR_LAGRANGE_MULTIPLIER, self.main_model_part.Nodes, label, 0)
Expand All @@ -510,24 +539,24 @@ def _get_if_is_interval(self):
self -- It signifies an instance of a class.
"""
current_time = self.main_model_part.ProcessInfo[KM.TIME]
if(self.interval.IsInInterval(current_time)):
if self.interval.IsInInterval(current_time):
return True
else:
return False

def _compute_number_nodes(self):
# We compute the number of nodes of the geometry
if (self.predefined_master_slave is True and self.dimension == 3):
if self.predefined_master_slave is True and self.dimension == 3:
slave_defined = False
master_defined = False
for cond in self.computing_model_part.Conditions:
if (cond.Is(KM.SLAVE)):
if cond.Is(KM.SLAVE):
number_nodes = len(cond.GetNodes())
slave_defined = True
if (cond.Is(KM.MASTER)):
if cond.Is(KM.MASTER):
number_nodes_master = len(cond.GetNodes())
master_defined = True
if (slave_defined and master_defined):
if slave_defined and master_defined:
break
else:
number_nodes = len(self.computing_model_part.Conditions[1].GetNodes())
Expand All @@ -543,17 +572,17 @@ def __get_integration_area(self):
"""
# We compute the total integrated area, for debugging
total_area = 0.0
if (self.dimension == 2):
if self.dimension == 2:
exact_integration = KM.ExactMortarIntegrationUtility2D2N(3)
else:
number_nodes, number_nodes_master = self._compute_number_nodes()
if (number_nodes == 3):
if (number_nodes_master == 3):
if number_nodes == 3:
if number_nodes_master == 3:
exact_integration = KM.ExactMortarIntegrationUtility3D3N(3)
else:
exact_integration = KM.ExactMortarIntegrationUtility3D3N4N(3)
else:
if (number_nodes_master == 3):
if number_nodes_master == 3:
exact_integration = KM.ExactMortarIntegrationUtility3D4N3N(3)
else:
exact_integration = KM.ExactMortarIntegrationUtility3D4N(3)
Expand Down Expand Up @@ -598,12 +627,12 @@ def __generate_search_model_part_from_input_list(self, param, key = "0"):

# We set the interface flag
KM.VariableUtils().SetFlag(KM.INTERFACE, True, partial_model_part.Nodes)
if (len(partial_model_part.Conditions) == 0):
if len(partial_model_part.Conditions) == 0:
KM.Logger.PrintInfo("Contact Process", "Using nodes for interface. We recommend to use conditions instead")
else:
KM.VariableUtils().SetFlag(KM.INTERFACE, True, partial_model_part.Conditions)

if (self.preprocess is True):
if self.preprocess is True:
id_prop = self.settings["search_property_ids"][key].GetInt()
if (id_prop != 0):
sub_search_model_part.SetProperties(self.main_model_part.GetProperties(id_prop))
Expand Down Expand Up @@ -635,7 +664,7 @@ def __detect_skin(self, model_part, key = "0"):
sub_search_model_part_name = "ContactSub"+key
self._get_process_model_part().CreateSubModelPart(sub_search_model_part_name)
detect_skin_parameters["name_auxiliar_model_part"].SetString(sub_search_model_part_name)
if (self.dimension == 2):
if self.dimension == 2:
detect_skin = KM.SkinDetectionProcess2D(model_part, detect_skin_parameters)
else:
detect_skin = KM.SkinDetectionProcess3D(model_part, detect_skin_parameters)
Expand All @@ -652,7 +681,7 @@ def __assume_master_slave(self, key = "0"):
key -- The key to identify the current pair
"""
# When all conditions are simultaneously master and slave
if (self.settings["assume_master_slave"][key].size() > 0):
if self.settings["assume_master_slave"][key].size() > 0:
self.predefined_master_slave = True
else:
self.predefined_master_slave = False
Expand Down